Are neutron stars crushed? Gravitomagnetic tidal fields as a mechanism for 

binary-induced collapse 



Marc Favatfjll 

Department of Astronomy, Cornell University, Ithaca, NY 14853 
(Dated: Received 23 October 2005; published 5 May 2006 in Phys. Rev. D 73 104005) 

Numerical simulations of binary neutron stars by Wilson, Mathews, and Marronetti indicated 
that neutron stars that are stable in isolation can be made to collapse to black holes when placed in 
a binary. This claim was surprising as it ran counter to the Newtonian expectation that a neutron 
star in a binary should be more stable, not less. After correcting an error found by Flanagan, 
Wilson and Mathews found that the compression of the neutron stars was significantly reduced 
but not eliminated. This has motivated us to ask the following general question: Under what 
circumstances can general-relativistic tidal interactions cause an otherwise stable neutron star to 
be compressed? We have found that if a nonrotating neutron star possesses a current-quadrupole 
moment, interactions with a gravitomagnetic tidal field can lead to a compressive force on the star. 
If this current quadrupole is induced by the gravitomagnetic tidal field, it is related to the tidal 
field by an equation-of-state-dependent constant called the gravitomagnetic Love number. This is 
analogous to the Newtonian Love number that relates the strength of a Newtonian tidal field to the 
induced mass quadrupole moment of a star. The compressive force is almost never larger than the 
Newtonian tidal interaction that stabilizes the neutron star against collapse. In the case in which a 
current quadrupole is already present in the star (perhaps as an artifact of a numerical simulation) , 
the compressive force can exceed the stabilizing one, leading to a net increase in the central density 
of the star. This increase is small (< 1%) but could, in principle, cause gravitational collapse in a 
star that is close to its maximum mass. This paper also reviews the history of the Wilson-Mathews- 
Marronetti controversy and, in an appendix, extends the discussion of tidally-induced changes in 
the central density to rotating stars. 
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I. INTRODUCTION AND SUMMARY 

Binary systems of two neutron stars (NSs) or a neu- 
tron star and a stellar-mass black hole (BH) are possible 
sources of gravitational waves (GWs) for current Q, 
and future |3( GW interferometers. To extract informa- 
tion from these waves the stages of the coalescence must 
be modelled accurately. When the binary separation d is 
large (such that d^> R, where R is the radius of the NS), 
analytic post-Newtonian (PN) methods |j| can describe 
the binary dynamics accurately enough to allow detection 
and parameter extraction. However, as the binary sepa- 
ration decreases, the PN approximation (which assumes 
weak gravity and slow motion) becomes less and less ac- 
curate. At some point the system must be modelled by 
numerical simulations that account for strong gravita- 
tional fields and hydrodynamic effects. Several groups 
have developed numerical codes to simulate NS/NS sys- 
tems (e.g., see @ and Refs. 7-15 of :6{]). The detec- 
tion of GWs from NS/NS coalescences could yield in- 
formation about the equation of state (EOS) of ultra- 
dense nuclear matter, and about short-duration gamma- 
ray bursts 0, HI • Accurate predictions of the GW signal 
will be important for these purposes. 

Wilson, Mathews, and Marronetti (WMM) [![l(j were 
one of the first groups to simulate the hydrodynamics of 
NS/NS mergers in general relativity. Their simulations 
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made the surprising prediction that relativistic effects 
can compress neutron stars that are near their maximum 
mass, initiating collapse to black holes prior to the onset 
of the dynamical orbital instability that causes the stars 
to plunge and merge. This prediction, which is referred to 
by some as "star-crushing" or "binary-induced collapse," 
was highly controversial and ran counter to intuition ob- 
tained from the Newtonian result that a NS in a binary 
is more stable against collapse [Tlj . If true, this collapse 
instability would have important implications for the de- 
tection of NS/NS binaries using matched filtering. The 
energy loss from the collapse process would change the 
orbital phase and introduce additional EOS-dependent 
parameters in the inspiral waveform templates. Over 15 
papers appeared in the literature refuting WMM's claim. 
Details of this controversy are reviewed in Sec. ITT1 below 
and in Q ■ Kennefick |l2j provides a very interesting and 
readable account of the controversy from a sociological 
viewpoint. The WMM controversy largely subsided once 
Flanagan ^| discovered an error in one of WMM's equa- 
tions. Although correcting this error caused a substantial 
decrease in the crushing effect, some compression of the 
neutron stars remained 

Various an^ticlll[i]|^ Hj and nu- 

merical m m iioiiiirsinia i3i studies h ave 

claimed to rule out the star-crushing effect. However, 
none of these studies considered certain post-Newtonian, 
velocity-dependent tidal couplings or they constrained 
the NS velocity field to be either initially vanishing, coro- 
tating (where the NSs are rigidly rotating at the orbital 
frequency), irrotational (the NS fluid velocity has van- 
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ishing curl) , or described by ellipsoidal models (in which 
the velocity field is a linear function of the distance from 
the star's center of mass); see Sec. |H] for further discus- 
sion. These approximations have left open a loophole 
in the demonstration that the central density of a neu- 
tron star should always decrease when placed in a binary 
system. Specifically, there remains the possibility that 
gravitomagnetic tidal interactions could couple to com- 
plex velocity patterns inside a neutron star, causing the 
central density to increase. The purpose of this paper is 
to investigate whether such a mechanism can explain the 
residual compression observed in WMM's revised simu- 
lations 14] and, more importantly, to address the follow- 
ing general question: Are there any circumstances under 
which general-relativistic tidal forces can compress a neu- 
tron star? 

We find that there is a compression effect which can be 
briefly summarized as follows: In addition to the familiar 
Newtonian tidal field of its companion, the fluid of each 
NS also interacts with a gravitomagnetic tidal field gener- 
ated by the motion of its companion. If the NS fluid has a 
nonzero current-quadrupole moment, velocity-dependent 
tidal forces can lead to compression of the star, increasing 
its central density in certain circumstances and making 
it more susceptible to gravitational collapse. 

To describe this mechanism in mathematical language, 
begin by considering a nonrotating neutron star with 
mass M and radius R interacting with the tidal field of a 
binary companion with mass M' a distance d away. In- 
troduce the dimensionless book-keeping parameters e — 
M/R (which parameterizes the strength of the NS's in- 
ternal gravity) and a = R/d (which parameterizes the 
strength of tidal forces) . We use units with G = c = 1 . 
For our purposes, we can treat the star's internal self- 
gravity as Newtonian (see Appendix 0. Then at lead- 
ing order in e and a, the metric in the vicinity of the star 
with mass M can be expanded as 

.goo = -1 - 2$ - 2$ cxt + 0(e 2 ) + Q{ea 4 ) , (1.1a) 



.go, = C xt + O(6 3 / 2 a / 2 ) + 0{^a 7 ' 2 ) , (1.1b) 

9ij = (1 - 2$ - 2$ cxt )<% + 0(e 2 ) + 0(ea 4 ) , (1.1c) 

where $ = 0(e) is the star's self-gravitational Newtonian 
potential, and $ cxt = 0(ea 3 ) and C, cxt = 0(e 3/2 a 7/2 ) are 
the Newtonian and gravitomagnetic potentials describing 
the external tidal field. Inside and near the star these 
potentials satisfy a subset of the first post-Newtonian 
(1PN) Einstein field equations, V 2 $ ext = V 2 Cf 3rt = 
and V 2 <& = 47T/9, where p is the NS's mass density. Our 
metric expansion l|l.lfl is not a complete 1PN expansion 
but only includes Newtonian and gravitomagnetic terms. 
A detailed justification of the expansion l|l.l|l is given in 
Appendix None of the terms that we neglect affect 
our final results. 



The external potentials in (|1.1H can be expanded as 
power series in the spatial coordinates x % whose origin 
follows the star's center of mass worldlinc: 

$ oyt = ^£abX a x b + 0(x 3 ) , (1.2) 

= -\e im W> lX ^x l +0{x 3 ) , (1.3) 

where £ij(t) and Bij(t) are electric- type and magnetic- 
type tidal moments. These moments are symmetric 
and trace-free (STF) tensors. They can be written in 
terms of the Riemann tensor of the external (tidal) pieces 
of the metric (|1.1J1 evaluated at the spatial origin via 
£ij(t) = Roioj and Bij(t) = \ei pq R pq jo- See Appendix 1X1 
for further discussion. 

In addition to the Newtonian tidal force, magnetic- 
type tidal fields introduce acceleration terms in the hy- 
drodynamic equations that resemble the vector-potential 
and Lorentz-force terms from electromagnetism, 

a cxt = -V$ oxt - C° xt + v x B . (1.4) 

Here B — V x £ cxt is the gravitomagnetic field, v is 
the internal fluid velocity measured with respect to an 
inertial frame who's origin coincides with the star's center 
of mass, and an overdot denotes a time derivative. 1 As we 
will show below ( Sees. ITTT1 and IIV|) . gravitomagnetic tidal 
forces can compress a star if the angle average of the vxB 
Lorentz-like force is nonzero and inward pointing. Such 
a force can only arise if the star's internal velocity field 
has a component in the subspace spanned by the I = 2 
magneticlike vector spherical harmonics Y B,lm oc x x 
Vy !m . (See Appendix iBl or Thorne 29] for a discussion 
of vector spherical harmonics.) The velocity field will 
have a nonzero component of this type if and only if the 
star's current-quadrupole moment Sij is nonzero. In the 
weak-field, slow-motion limit the current quadrupole is 
defined by 

Sij = J X(i€j) ab x a pVb d 3 x , (1.5) 

where p is the mass density, is the fluid velocity, and 
the parentheses denote symmetrization. Such a velocity 
field is depicted in Figures and [21 

If the gravitomagnetic tidal field Bij it) is slowly vary- 
ing, and if the star is initially static, then the £ ext term 
in Eq. i|1.4|) induces a velocity field given by v = — £ oxt . 
The corresponding current-quadrupole moment is 

Sij = T*MB*Bij . (1.6) 

Here 72 is the gravitomagnetic Love number, a dimension- 
less constant that depends on the NS equation of state 



1 We have dropped other 1PN terms from the tidal acceleration. 
This is justified in Appendix [X] Retaining them does not affect 
our results. 
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FIG. 1: (color online) Internal velocity field v of a nonrotating 
neutron star with a current-quadrupole moment induced by 
the tidal field of an orbiting companion. The arrows denote 
the velocity vectors and are generated by Eqs. 13. 1H and 
l|3.35|l . The velocity field is stationary in a coordinate frame 
rotating at the binary's orbital angular velocity (which points 
perpendicular to the equatorial plane of the star). 



(see Sec. HUB and Appendix EJ. 2 This process is analo- 
gous to the Newtonian tidal distortion of stars, wherein 
the electric-type tidal field Eij induces a mass quadrupole 
moment 2^ given by 

lij = -\k 2 R 5 E %3 . (1.7) 

Here k 2 is the dimensionless Newtonian Love number (see 
chapter 4.9 of [31| or Appendix IT)1). 

As shown in Sec. II I II the gravitomagnetically induced 
velocity field (Figures and |2J drives the fundamental 
radial mode of the NS (along which compression and 
decompression occur) via a combination of the Lorentz 
v x B and nonlinear advection (v ■ V)i> terms. (Figure|3| 
shows the total gravitomagnetic tidal acceleration acting 
on the fluid in an inertial reference frame whose origin in- 
stantaneously coincides with the NS center of mass.) Up 
to order 0(a 7 ), the resulting change in central density is 



2 This gravitomagnetically-induced current-quadrupole moment is 
related to Shapiro's gravitomagnetic induction of circulation 
in a NS by the gravitational field of a spinning black hole; see 
Sec. IHI Bl and Appcndix|C] The ellipsoidal model of a NS used in 
30] excluded current-quadrupole moments. Our analysis is also 
applicable to a spinning BH or any other source that produces a 
Bij tidal field. 




FIG. 2: (color online) Same as Fig. Q but showing only the 
induced velocity field on a slice through the equatorial plane. 
The velocity current loops set up in the star are easily seen. 



= Cl £y (*)£«(*) + c 2 Bij(t)B ij (t) , (1.8) 

Pc 

where the constants c\ and c 2 have units of [length] and 
depend on M, R, and the equation of state 01 . In a 
binary, the tidal fields scale as £y ~ M'/g? 3 and Bij ~ 
(M'/d 3 )y/(M + M')/d, so the two terms scale as 0(a 6 ) 
and 0(a 7 ), respectively. The first term in Eq. I|1.8f) is the 
Newtonian tidal-stabilization term. Its sign (c± < 0) has 
been computed for relativistic stars by Thorne |17| ; Lai 
[TT| and Taniguchi and Nakamura have computed 
its value for Newtonian stars, ci » — 0.38i? 6 /M 2 (for a 
r = 2 polytrope) . Its derivation is reviewed in Appendix 
El One of the main results of this paper is the mag- 
nitude and sign of the coefficient c 2 : it is positive and 
has the value c 2 « 0.064R 5 /M (also for a T = 2 poly- 
trope). This term therefore tends to compress the star. 
However, its size is not large enough to overcome the 
decompressive effect of the first term. Therefore, non- 
rotating neutron stars with no preexisting velocity fields 
suffer no net compression when placed in a binary. In 
Appendix [H] we briefly discuss how to extend our results 
to rotating stars. 

In Sec. lIVI we consider the possibility that the neutron 
star is not initially unperturbed but instead has a pre- 
existing current quadrupole. (By "preexisting" we mean 
that the current quadrupole does not arise through the 
mechanism of gravitomagnetic tidal induction discussed 
here.) Viscosity will damp astrophysical sources of a cur- 
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FIG. 3: (color online) The total gravitomagnetic tidal ac- 
celeration of a nonrotating neutron star interacting with a 
binary companion. Only a slice through the equatorial plane 
is shown. The arrows represent the acceleration vectors in an 
inertial reference frame centered on the star's center of mass 
and are given by Eqs. 13.181 and 13.351 . Their magnitude 
vanishes at the center of the star. As in Figs. 0and|5|the ac- 
celeration field rotates at the binary's orbital angular velocity. 
The radially-inward pointing acceleration indicates that the 
gravitomagnetic tidal interaction causes compression. 



rent quadrupole on a timescale t v ; s < 1 day. 3 Any ve- 
locity currents arising from the formation of the NS will 
be damped long before the NS comes close to merging. 
Unless they are generated shortly before coalescence, as- 
trophysical preexisting current quadrupoles are unlikely. 
However, a current quadrupole could be present as a nu- 
merical artifact in a NS/NS simulation. Approximations 
to the equations of motion, numerical errors, artificial vis- 
cosity, or the method of choosing the initial data could 
possibly lead to a nonzero current-quadrupole moment. 
It is possible that such a numerical artifact was present 
in the WMM simulations 01 ■ I n an y case, the presence 
of a preexisting current quadrupole affects the change in 
central density by replacing the second term in Eq. (|1.8|) 
with 45y(t)B«(t), where c' 2 » 2.9R/M 2 (for a T = 2 
polytrope). This term scales like a 7 / 2 , and at large sep- 
arations it actually dominates over the Newtonian tidal- 
stabilization term. The time dependence and sign of this 



3 The viscous time is T v j a ~ pR 2 /r), where p is the density and 
n ~ 4.7 x 10 19 g cm _1 s _1 (T/10 8 K) _2 |>/(2.8 x 10 14 g cm" 3 )] 2 
is the coefficient of shear viscosity (this is valid at low tem- 
peratures T < 10 9 when protons and neutr ons are superfluid 
and electron-electron scattering dominates |32| ). This gives 
T vis ~ 0.0019 day(R/10km) 2 (T/10 6 K) 2 (p/10 15 g cm" 3 )" 1 



term depends on the unknown functional form of Stj (t) . 
If we assume that Sij(t) is constant, the 0(a 7 / 2 ) term 
oscillates in sign at the orbital period, and a net com- 
pressive force results during parts of the orbital phase. 
For plausible values of Sij(t), the net change in central 
density is small for Newtonian stars, < 1% (see Fig. [JjJ, 
but it could be large enough to cause collapse if the NS 
is close to its maximum mass. 

In the remainder of this article and in its appendices, 
we provide the details of the analysis summarized above. 
But first we give further motivation for our analysis by 
reviewing the history of the WMM star-crushing contro- 
versy (Sec. llTjl . 

Throughout this paper we follow the notations and 
conventions of Misner, Thorne, and Wheeler [sli ] 
(MTW). We assume geometric units with G = c = 1. 
Time and space coordinates are denoted by x a = (t,x J ). 
Spatial indices (in a Cartesian basis) are raised and low- 
ered using dij. Repeated spatial indices are summed, 
whether or not they are up or down. Spatial partial 
derivatives are denoted by V,; and time derivatives are 
denoted by an overdot, / = df/dt. Spacetime indices 
and covariant derivatives are rarely used. 

II. A BRIEF HISTORY OF A CONTROVERSY 

To help motivate our analysis, it is useful to review 
the history of the star-crushing controversy, focusing on 
the arguments for and against crushing and the approxi- 
mations that are assumed in the various arguments. We 
begin by reviewing the ori gina l WMM simulations. 

WMM's simulations |gL Il0| relied on two important 
assumptions, which have come to be called the Wilson- 
Mathews approximation 4 : First, the spatial metric sat- 
isfies the spatial conformal flatness (SCF) condition, 
Jij = 4 > Sij, where jij is the 3-metric of a spacelike hy- 
persurface and <f> is the conformal factor. The SCF con- 
dition simplifies the form of the hydrodynamic and field 
equations, neglects gravitational radiation in the spatial 
3-metric, and is generally accurate only to 1PN order. 
However, it is exact for situations with spherical sym- 
metry and very accurate for rapidly rotating relativistic 
stars |3Gll . Although widely used by many groups, the 
SCF condition was suspected by some to be the source 
of WMM's crushing effect (but see Sec. HEl below). The 
second assumption is a quasiequilibrium approximation 
in which the terms involving the time derivatives of the 
gravitational degrees of freedom (the spatial metric 7^ 
and extrinsic curvature ify) are dropped from the equa- 
tions of motion. This is thought to be a good approx- 
imation at large separations when GWs hardly modify 
the orbital dynamics. Combined with the SCF condi- 
tion, this assumption reduces the equations for the grav- 



4 A self-contained descriptio n of the WMM simulations is also 
found in their rec ent book |34|| . For a shorter review of then- 
work, see Ref. l35l . 
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itational field to flat-space elliptic equations. Given an 
initial matter distribution, WMM first solve the momen- 
tum and Hamiltonian constraint equations for the grav- 
itational field. The hydrodynamics equations (coupled 
to the gravitational field) are then evolved to the next 
time slice. Instead of also evolving the gravitational field 
variables, the constraint equations are solved again at 
that time slice and the process is iterated. Gravitational 
waves are calculated via a multipole expansion and their 
effect on the neutron stars is accounted for by adding a 
radiation-reaction potential to the hydrodynamics equa- 
tions. WMM also employ what they refer to as a "re- 
alistic equation of state" . This zero-temperature, zero- 
neutrino-potential EOS |l0ll37t l38l| is softer (smaller val- 
ues of M/R) than polytropic equations of state used by 
other groups and shows greater compression in their sim- 
ulations. This EOS was motivated by matching models 
of SN 1987A to the observed neutrino signal [39 ]. 

Unlike most other NS/NS simulations of that time, 
the WMM simulations used unconstrained hydrodynam- 
ics — they did not constrain the binary to be corotating 
or irrotational. Even though more recent NS/NS sim- 
ulations also use unconstrained hydrodynamics (JjJ and 
references therein), they all constrain the stars in their 
initial data sets to be either corotating or irrotational. 
In the WMM simulations, the initial data is formulated 
differently (Sec. Ill of 10]): An initial " guess" solution 
from the Tolmen-Oppenheimer-Volkoff equation for each 
star is placed on the grid in a corotating configuration. 
The stars are then allowed to relax to an equilibrium 
configuration. This is accomplished by solving the field 
equations and evolving the hydrodynamics without radi- 
ation reaction. An artificial damping of the fluid motion 
is imposed and slowly removed as stars reach an equilib- 
rium state. The resulting equilibrium configurations are 
neither corotational nor irrotational, but result in stars 
with almost no intrinsic spin (Sec. IV E of 01). As dis- 
cussed in Sec. lIVI this method of choosing the initial data 
sets could possibly be the source of compression. 

The main result of the initial WMM simulations was 
that initially stable neutron stars could be highly com- 
pressed: a star ~ 12% from its maximum mass has a 
change 5p c of its central density p c given by 5p c /p c w 0.51 
at a proper separation of 68 km [14(. The simulations 
indicated that the central density increased according 
to p c oc U A , where U 2 — UiU 1 and Ui are the spatial 
components of the 4- velocity (see Fig. 2 of HJj). WMM 
also found that the binary's orbit would become unstable 
at an orbital separation that was larger (by a factor of 
~ 1.4) than the PN prediction 0. 

In the years following WMM's initial publications, sev- 
eral papers appeared claiming that neutron stars in a 
binary should be stabilized and not compressed. These 
were followed by a rebuttal paper by WMM Lai 
|ll| used an energy variational principle (including 1PN 
corrections to the star's self-gravity) to show that a New- 
tonian tidal field decreases the central density according 
to Sp c /p c = -2.7(M'/M) 2 (R/d) 6 (for an n = 3/2poly- 
trope at its maximum mass in isolation; see also [l(j and 
Appendix [E] of this paper) . Wiseman showed that 



there was no change in central density in a binary at 1PN 
order, but he neglected tidal effects. Brady and Hughes 
[TH examined a point particle with mass p <C M orbiting 
a static, spherical NS and showed that there is no change 
in central density at linear order in p. Thorne [l7| showed 
that fully-relativistic, static or rotating NSs are stabilized 
by an electric-type tidal field. Although these papers 
(lBL UtI fl8| consistently applied their approximations, 
they did not include the velocity-dependent forces that 
WMM attribute their compression to |4J|, an d they did 
not consider the gravitomagnetic interactions that we in- 
vestigate here. Shibata and Taniguchi |4l| and Lombardi 
et al. floj both examined equilibrium sequences of com- 
pressible ellipsoids at 1PN order. Shibata and Taniguchi 
considered corotating binaries while Lombardi et al. con- 
sidered corotating and irrotational ones. Both also found 
that the NSs were stabilized, but WMM claim that they 
also ignored the relevant velocity-dependent terms |4C| . 
Shibata et al. j2|| performed 1PN hydrodynamics simula- 
tions for corotating and irrotational binaries and also saw 
no signs of compression. WMM speculated that this was 
due to the unrealistically soft EOS (with M/R w 0.023) 
used by Shibata et al. |2g. For the very close separations 
examined in that paper, WMM claimed that the tidal 
stabilization overwhelms any compression effect E9. 

In a series of papers, Baumgarte et al. 0, |22i 
simulated corotating NS/NS binaries using the SCF and 
quasiequilibrium conditions, finding that the stars were 
stabilized. However, their simulations did not contra- 
dict the WMM results since the centrifugal force tends to 
stabilize the star in corotating binaries. Further, WMM 
showed analyti cally that their compression effect vanishes 
for corotation 40]. This indicated to WMM that the 
compression was probably due to the nonrigidly-rotating 
motion set up in the NS fluid [40|. 

A matched-asymptotic-expansion analysis of the 
crushing effect was performed by Flanagan |16|. He 
showed that, to all orders in the strength of internal grav- 
ity of each NS, the leading-order terms in a tidal expan- 
sion of the change in central density are given by Eq. 11. 8|) 
above. Flanagan also showed that the coefficient c\ of the 

leading 0(a 6 ) term has the form F [1 + F x e + F 2 e 2 H ] , 

where Fq, F\, F2, . . . are constants. His analysis did not 
determine the overall sign of c\, but he concluded that 
since Fq was shown by Lai's [ll| Newtonian analysis to be 
negative, the central density of a NS in a binary will de- 
crease unless i<2, • • • are negative and large. Thome's 
\Ti\ relativistic analysis showed that the entire coefficient 
ci is negative, thus excluding the possibility of a sign 
flip. Flanagan did not determine the sign or magnitude 
of the coefficient C2 in Eq. i|1.8|) , which is one of the main 
results of this paper (although, in contrast to Flanagan, 
the internal gravity of each NS is Newtonian in our treat- 
ment). Flanagan's analysis accounts for gravitomagnetic 
tidal fields and velocity-dependent corrections to the hy- 
drodynamics that are induced by tidal interactions. It 
neglects, however, any crushing that could be caused by 
preexisting velocity fields. WMM indicate that such ve- 
locity fields may be responsible for their observed com- 
pression (43. We address this in Sec. IIVI 
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Despite the numerous claims that NSs in binaries are 
stabilized against collapse, there are a few analyses that 
hint that the binary-induced collapse of compact objects 
is possible. Shapiro ji^ considered a system of a "com- 
pact object" made up of a test particle in a close orbit 
around a nonrotating BH, perturbed by the Newtonian 
tidal field of a distant binary companion. Although the 
test particle has a stable orbit in isolation, the tidal field 
could cause the test particle to plunge into the BH. Duez 
et al. m| extended this analysis to a swarm of parti- 
cles and included relativistic effects neglected by Shapiro, 
confirming his conclusions. Alvi and Liu |44j also exam- 
ined the stability of a swarm of test particles but included 
the effects of magnetic-type tidal fields. They found that 
including magnetic-type tidal fields did not strongly af- 
fect the average radius of the cluster, but it did destabi- 
lize individual particles that were stable in the absence 
of magnetic-type tidal fields. 

Despite indications of binary-induced collapse, it seems 
unlikely that these models are relevant to situations 
where hydrodynamic forces are present. For circular 
orbits, the test particles in these simulations lie at 
the stable minimum of the effective potential of the 
Schwarzschild geometry (see chapter 25 of MTW, espe- 
cially Fig. 25.2). For particles close to the last stable 
orbit, this minimum is only marginally stable. The ex- 
ternal tidal forces perturb the test particles about this 
minimum. The direction and size of the perturbing tidal 
force depends on the relative orientation and separation 
of the particle and the tidal field. When the tidal per- 
turbation is small the particle rolls "up the hill" of the 
potential and then rolls back to the stable minimum. But 
if the tidal perturbation is large enough, the particle can 
be forced over the local maximum of the potential, caus- 
ing it to plunge into the BH's event horizon. Adding 
additional (magneticlike) tidal fields simply provides an 
additional force that will cause more particles to become 
unstable. The binary-induced collapse of a star is differ- 
ent because pressure and not orbital angular momentum 
supports the star against collapse. Collapse can only oc- 
cur if the angle average of the tidal force points radially 
inward. This is harder to achieve than accelerating a 
single particle to smaller radii. 

The controversy appeared to be resolved when Flana- 
gan 0] found an error in one of WMM's equations and 
showed that this error could account for the observed 
compression. The error was an incorrect definition of the 
momentum density in the momentum constraint equa- 
tion. Wilson and Mathews [1J| corrected this error and 
showed that the compression was reduced (by about a 
factor ~ 10) but not eliminated. (They also noted that 
the frequency of the last stable orbit moved closer to the 
post-Newtonian value.) For a T = 2 polytrope, Sp c /p c 
was reduced from 0.14 (at a 138 km separation) to 0.008 
(at a 118km separation). Using their realistic EOS and 
stars with a gravitational mass of 1.39Mq, Sp c /p c was 
reduced from 0.51 (at a 68 km separation) to 0.03 (at 
a 61km separation). Stars closer to the last stable or- 
bit showed a compression of ~ 10% but did not collapse 
(as they did in the uncorrected simulations; see Figure 



0}. But for stars close to their maximum mass (< 9% for 
their realistic EOS) and for very close (but stable) orbits, 
collapse to BHs could still occur. (In this case the coor- 
dinate separation between the stars was 2.4 times their 
coordinate radii.) The p c oc U 4 scaling of the central 
density also remained in their revised simulations. Wil- 
son and Mathews also state that the question remains 
as to whether their residual compression "is real or an 
artifact of the numerics" or the SCF approximation |14j . 

Wilson and Mathews continue to identify the observed 
compression as arising from enhanced self-gravity terms 
proportional to the square of the fluid velocity 0j 03 • 
These terms originate from the T I1 ^\T^ X term of the hy- 
drodynamics equations (here is a connection coef- 
ficient and T^ A is the energy- momentum tensor). Al- 
though they claim that tidal effects do not cause the 
compression ^tji the conventional understanding of the 
equivalence principle suggests that all gravitational in- 
teractions of a NS with an external body are tidal inter- 
actions. The analyses of Thorne 0J an d Flanagan 0] 
support this argument, as does the present paper. This 
suggests that the residual compression in [l4| might be an 




d [km] 

FIG. 4: (color online) Change in central density for the revised 
Wilson-Mathews simulations as a function of proper distance 
between the neutron star centers. The lower set of points 
(Table III of Q) corresponds to stars that are 12% from 
their maximum mass (in isolation). These stars were com- 
pressed, but did not collapse before the last stable orbit. A 
fit to these points (dashed, blue curve) indicates the scaling 
Spc/pc oc a 2 ' . The upper set of points (Table IV of Q ) cor- 
responds to stars that are 8.6% from their maximum mass (in 
isolation). In this case the stars did collapse (not shown here) 
and the orbit remained stable. A fit to these points (solid, 
red curve) indicates the scaling 5p c /p c oc a 1 ' 4 . In both cases, 
a "realistic" equation of state was used (see [Hill). 
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artifact of the computational scheme they have chosen. 
If the revised Wilson-Mathews [TJJ simulations contain 
some fluid circulation in their initial data, compression 
could occur via the mechanism discussed in Sec. IIVI be- 
low. We also note that, despite skepticism of their com- 
pression effect, Wilson and Mathews continue to invoke 
it as a mechanism to explain gamma-ray bursts [39t |46| 
and, recently, to propose a new class of Type I supernovae 



A. Compression in irrotational simulations 

Because of the controversial nature of their results, 
WMM developed an independent numerical code using 
the irrotational approximation |47|. (The hydrodynam- 
ics was unconstrained in their previous simulations.) In 
the irrotational approximation the fluid vorticity is zero. 
These simulations also show a small increase in central 
density (1.5% at 30 km separation for a T — 2 polvtrope) 
that is larger than the numerical errors estimated in |4jj 
but is within the possible error induced by the SCF con- 
dition. WMM 47] also claim that this compression is 
consistent with the irrotational simulations of Bonazzola 
et al. 0>E3- I n this section, we review the results of irro- 
tational NS/NS simulations from two independent groups 
which show no evidence for compression. This indicates 
that the small compression seen in WMM's irrotational 
simulations is unphysical. The observed compression is 
possibly due to the inaccurate treatment of a boundary 
condition or insufficient grid resolution. 

Other numerical groups have shown that no central 
compression occurs for NS /NS binaries in the irrotational 
approximation. Although WMM claim that Bonazzola et 
al. also see a small compression of order < 0.3% 

(see Figs. 12 and 13 of 0]), the central density decreases 
with decreasing orbital separation in their simulations (in 
contrast to WMM 47]) and is within the error induced 
by the SCF approximation. Uryu et al. [53, H3] also per- 
formed irrotational simulations and see a decrease in cen- 
tral density at small separation. While they also see oscil- 
lations in which the central density increases by ~ 0.5% 
(see Fig. 6 of 28]), they claim that this is due to the 
errors of their finite-difference scheme and of their Leg- 
endre expansion of the gravitational field (Sec. Ill D of 
|28j). Furthermore, after improving their method of de- 
termining the stellar surface, the slight increase in central 
density seen in 0, [13 is removed and the central den- 
sity decreases monotonically (by ~ 1%) with decreasing 
separation (see Fig. 2 and footnote 3 of Taniguchi and 
Gourgoulhon |U; see also Figs. 12-14 of Hj). 



More precisely, the specific momentum density per baryon is ex- 
pressed as the gradient of a potential, hu^ = VuW, where 
is the 4-velocity, V u i s a covariant deri vati ve, and h is the rel- 
ativistic enthalpy l47l . See Teukolsky |4Sl for a discussion of 
the irrotational approximation in NS/NS simulations; see also 
Appendix IH1 of this paper. 



The source of compression in WMM's irrotational sim- 
ulations 1S m °st likely not the SCF o r q uasiequi- 
librium assumptions. The French |H Hfl H3 and 
Japanese [53, numerical groups also make these as- 
sumptions but do not see compression. Further, Wilson 
[52| examined the head-on collision of two NSs using two 
separate simulations: one in full general relativity and 
the other using the SCF condition. He found similar 
levels of compression in both cases, indicating that the 
SCF condition is not a likely culprit. See Appendix B of 
Baumgarte and Shapiro 5] for a further discussion of the 
validity of the SCF condition. 

The primary difference between the irrotational sim- 
ulations of WMM and those of the other groups is the 
numerical technique used: The Japanese group used a 
multidomain, finite-difference method with surface-fitted 
spherical coordinates (which allow accurate resolution of 
the stellar fluid and surface). The French group used 
an even more accurate multidomain spectral method, 
also with surface-fitted spherical coordinates. WMM's 
technique is the least accurate: a single-domain finite- 
difference method with Cartesian coordinates. Both the 
French and Japanese groups point out a likely source of 
error in the WMM |47| simulations: the use of Cartesian 
coordinates and an approximate treatment of the bound- 
ary condition for the velocity potential \& that treats the 
stellar surface as spherical [see Eq. (19) of [5j| and the 
discussion in Sec. V A of [23 and Sec. VII A of |51|]. 
This issue is also discussed in Sec. 9.3 of 

It is also possible that low grid resolution is the source 
of the WMM compression [54( . Since they were using the 
best grid size possible at the time, it was not possible to 
estimate the error due to poor resolution in ^7|- Re- 
gardless of the precise source of error, the fact that more 
accurate simulations do not observe compression strongly 
suggests that the compression seen in WMM's irrota- 
tional simulations is unphysical. If poor grid resolution 
is the source of the compression in their irrotational sim- 
ulations, then it seems plausible that low resolution may 
also be the source of compression in the revised Wilson- 
Mathews simulations |14| using unconstrained hydrody- 
namics. However, we will also discuss in Sec. IIVI the 
possibility that the compression in [14( is related to the 
fact that the initial data sets in those simulations were 
neither corotational nor irrotational. 

Many numerical groups use the irrotational approxi- 
mation to cither simplify the evolution equations, or to 
determine the initial data when solving the constraint 
equations. The irrotational approximation is frequently 
motivated by the findings of Kochanek [55| and Bildsten 
and Cutler 5(| that the NS viscosities are too small to 
allow binaries to be tidally locked. However, this is more 
an argument against corotation than it is one in favor of 
irrotation. The irrotational assumption is also motivated 
by Kelvin's circulation theorem — in the absence of viscos- 
ity, initially irrotational flows remain irrotational; see Ap- 
pendix[U]for discussion. Irrotation is widely adopted pri- 
marily because it simplifies the hydrodynamic equations. 
However, there are physically well-motivated reasons to 
consider more general fluid configurations. Although re- 
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alistic NSs will not be corotating, they will have some 
intrinsic spin, thus violating the irrotation assumption. 6 
The much studied r-modes in rotating stars are another 
example of a velocity configuration that does not fit into 
the corotation or irrotation class. The excitation of these 
r-modes could lead to small effects on the GW signal, 
even in the low frequency (10Hz < / < 100Hz) regime 
|58[ . Although recent NS/NS simulations use uncon- 
strained hydrodynamics and full general relativity (see 
and references therein), they constrain the initial data 
sets to be corotational or irrotational. The WMM simu- 
lations |9l. ll0l.lT4| do not make this assumption. This pro- 
vides further motivation for our examination in Sec. lIVI of 
the coupling of preexisting current quadrupoles to tidal 
fields. 



III. GRAVITOMAGNETIC CONTRIBUTION 
TO THE CHANGE IN CENTRAL DENSITY 

A. Equations of motion 

To determine if a NS interacting with external tidal 
fields is compressed, we will compute the change in cen- 
tral density of the star by solving the fluid equations of 
motion. Begin by considering a star with mass M and 
radius R (in isolation) interacting with the external grav- 
itational field of a binary companion (characterized by a 
mass M' at a distance d). Assume that the star is ini- 
tially static in the following sense: when the binary sep- 
aration is very large the stellar fluid configuration is that 
of an unperturbed, nonrotating star in hydrostatic equi- 
librium. If one expands the metric in the local proper ref- 
erence frame of the star [as in Eqs. (|l.l|l ] and substitutes 
into the conservation of energy-momentum equation for 
a perfect fluid, the leading-order response of the star to 
the external gravitational field can be described by 

^+V i (pv i ) = 0, (3.1a) 



flit- V7 P 

^ + ( W fe V fe K = -^--V l $ + ar t , (3.1b) 



V 2 $ = 4irp. (3.1c) 

These are just the continuity, Euler, and Poisson equa- 
tions for a star with baryon density p, internal velocity Vi, 
pressure P, and Newtonian self-gravity $, augmented by 
an external driving force which is the 1PN point-particle 
acceleration (see chapter 9 of Weinberg j5^), 

-v 2 Vi$ ext + 4v 2 (v k V k )<i> cxt , (3.2) 



where B = Vx £ oxt . We also assume a barotropic EOS 
P = P{p). In the above equations we have ignored all 
PN corrections to the fluid equations except for the terms 
in the external acceleration a| xt . This is justified in Ap- 
pendix^] None of the terms that we drop will affect the 
leading-order corrections to the change in central density. 
The potentials <& cxt and C° xt that appear in Eqs. i|l.l|) 
and lE21) 

can be expressed in terms of the electric and 
magnetic- type tidal moments as in Eqs. 1|1.2[> and (|1.3fl . 
We will ignore the tidal octupole moment contribution, 
$ cxt = \£ijkX l x : > x k and higher moments; they will affect 
the central density at order 0(a 8 ) and higher. 

For our purposes we will only need to consider the 
lowest order tidal expansion of the first three terms in 
Eq. J3~2l : 

af xt = -eeSijxi + e B (^e ijk Bjx k x l - 2e ljk B kl v^ x l ^j 

+ 0(x 3 ), (3.3) 

where eg and eg are dimensionless book-keeping con- 
stants proportional to their respective tidal moments. 
They will be set to unity at the end of the calculation. 
One can explicitly show that for nonrotating stars, the 
terms in Eq. I|3.2[) that we have neglected will affect nei- 
ther the change in central density up to order 0(a 7 ) nor 
the leading-order contribution to the induced current- 
quadrupole moment; see Appendix lAl 

B. Second-order Eulerian perturbation theory 

To determine the influence of the external tidal fields 
on the structure of our star, we treat the tidal accelera- 
tion as a small perturbation whose size is parameterized 
by a dimensionless book-keeping parameter e. The den- 
sity, pressure, internal gravitational potential, and stellar 
velocity field are then expanded as 

p(t,x) = p(°)+ £ p(i)+ e y 2 )+... , (3.4a) 

P(t,x) = p(°> + eP w +e 2 P (2) +••• , (3.4b) 

<f>(t,x) = $ (0) +e$ (1) +e 2 $ (2) + ••• , (3.4c) 

v(t,x) = w<°) +£t>« +eV 2 ) +•■• , (3.4d) 

and substituted into the fluid equations (|3.1() . Each equa- 
tion is then solved order by order in e. For an initially 
static star v^ ' = 0. 7 

In a general analysis one could pick e = a and use the 
full expression for af xt in Eq. Q3.2p . However, one would 
find that the contribution to 5p c / p c at 0(a e ) would come 
solely from the leading-order Newtonian tidal term pro- 
portional to £ij, while the 0(a 7 ) contribution to Sp c /p c 



See Marronetti and Shapiro l57i for recent work that treats 
NS/NS binaries with arbitrary spin. 



For a slowly-rotating star in a tidal field, one would choose v^ ' = 
f2 X x and expand the fluid variables in both the tidal expansion 
parameter e and the angular velocity f2; see Appendix |H] for an 
analysis of this case. 
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and the leading contribution to Sij would only come from 
the gravitomagnetic terms in a| xt . It is therefore much 
simpler for our purposes to expand separately in either 
eg or eg. To compute the 0(a 6 ) tidal-stabilization term, 
one would set e = eg, es — in Eqs. I|3.1[l . i|3.3[l . an d 
(|3.4|) . and expand to 0{e\). This leading-order tidal- 
stabilization term is actually more difficult to compute 
than the 0(a ) destabilization term that we compute be- 
low. The tidal-stabilization term has also been computed 
by other methods [II], uM > this is reviewed in Appendix 
IE1 We will simply use the result of Taniguchi and Naka- 
mura [lj| for the change in central density of a T = 2 
Newtonian polytrope, 



5p c 

Pc 



45 fM'\ 2 /i? x 6 



2tt 2 V M 



(3.5) 



To compute the gravitomagnetic destabilization term, 
we set s = eg, e £ = in Eqs. (|3.ip . (|3.3|) . and l|3.4|l . 
expand, and solve the fluid equations order by order in 
eg. At order 0(e e ) we have the standard equations for 
a star in hydrostatic equilibrium, 



(0) 







(3.6) 



[Poisson's equation is also satisfied at each order in the 
expansion: V 2 $ ( ") =47rp( n ).] 
At order 0(sg) we have 



(3.7) 



and 



,(0) 



(i) 



0>f 



= -v i p( 1 )-p( 1 )v J $(°)-p(°)v l <i>( 1 )-p( )4 c 



(3.8) 

Combining the time derivative of (|3.7|) with the diver- 
gence of ((3.8j) and Poisson's equation gives 



gy*) 

<% 2 



8W + V 2 P« 

+V ! p (0) V l $ (1) + V lP (1) V l $ (0) , (3.9) 



where we have used V l [p^ %f xt ] = from Eq. p. 3(1 and 
= p( )(r). If our initial conditions state that there 
are no fluid perturbations at early times [so that at or- 
der O(eg) and higher, p, P, and v and their first 
time derivatives vanish as t — ► — oo], then the solution to 
Eq. O is 

(i) = p(i) = = o _ 



pvv = ^v^ = $vv =U . (3.10) 

Equation l|3.8|) then reduces to = — C° xt . In an inspi- 
ralling binary Bij — > as t — > — oo, and the leading-order 
velocity becomes 



= -cr = ^j k B^x k x i 



(3.11) 



This induced velocity shows that, in the absence of 
viscosity, a nonrotating star responds to the gravitomag- 
netic vector potential without resistance (like a spring 



with a vanishing spring constant). In rotating stars the 
Coriolis effect provides a restorin g for ce, and the gravito- 
magnetic field excites an r-mode |58|; the velocity (|3.11|) 
is the zero-rotation limit of the r-mode excitation. (See 
Figures ^ and |21 for a graphical depiction of this velocity 
field.) The velocity field l|3.11|l can be expressed as a sum 
of I — 2 magnetic-type vector spherical harmonics 



X 1 ) 



2 

E 

m=-2 



B 2 v m (.t,r)Y % 



B,2r. 



with 



B 2 v m {t,r) = - 



15 V 3 tiy% > 



(3.12) 



(3.13) 



(see Appendix 151 for definitions). Such a velocity field 
would be excluded by numerical simulations that en- 
force corotation. It would also be excluded in analy- 
ses that model each NS as an ellipsoid with an internal 
fluid velocity that is a linear function of the coordinates. 
However, such a velocity field would be permitted in a 
relativistic irrotational simulation. This seems puzzling 
at first because has nonvanishing Newtonian vortic- 
ity, LOi = [V x = 2Bijxi . The resolution is that 

the 1PN limit of the relativistic irrotational condition, 
V x (v^ + £ cxt ) — 0, is satisfied [3(J. In contrast, a ro- 
tating star or a nonzero frequency r-mode would not sat- 
isfy the relativistic irrotational condition. See Appendix 
El for further discussion. 

The velocity field (|3.11|) endows the NS with an 
induced current-quadrupole moment. Substituting 
Eq. H3.11fl into Eq. (|1.5|) gives Sy = jBij , where 
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8tt 
15 



pr 6 dr = j 2 MR 4 



(3.14) 



and 72 is the gravitomagnetic Love number. For a uni- 
form density Newtonian star, 72 = 2/35; for a T = 2 
polytrope 72 = 2(n 4 - 20n 2 + 120)/(15tt 4 ) ps 0.0274 (see 
Appendix . 

At order O(eg) we have the equations necessary to 
compute the change in central density at 0(a 7 ): 



djP_ 

dt 



(2)l 



0, 



(3.15) 



and 



,(0) 



Vi$< 2 > + 7(ojVi$(°) = a?* , (3.16) 



P 



dt py 
where 

.tot _ x B] . _ . v]v (l) 



(3.17) 



This acceleration term shows that the second-order 
perturbations are driven by a combination of the 
Lorentz-type gravitomagnetic and the nonlinear convec- 
tive derivative terms. (See Figure|3|for a graphical depic- 
tion of a* ot .) Both terms are generated by the first-order 
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velocity perturbation to the star. Using Eq. 1)3. 11[) . the 
acceleration term a* ot can be expressed explicitly in terms 
of the gravitomagnetic tidal field, 



In the absence of external driving (a* ot = 0) , employing 
the standard e~ luJo:t ansatz for the mode time dependence 
yields the eigenvalue equation for the modes, 



(3.18) 



where 
Hijki 



^ BikBjl — BakBalSij — -ecajtiblBakBbt 



(3.19) 



C. Radial Lagrangian perturbations 



To compute the change in central density, we first note 
that, since the first-order perturbations to the density, 
pressure and self-gravity vanish, Eqs. (|3.15() and (|3.16() 
can be recast as the equation for a linear Lagrangian 
perturbation of an initially static star. This is done by 
relabelling p = p<-°\ P = P(°\ $ = ¥°\ Sp = p& , 
SP = P (2) , 5$ = <f> (2) , and 5v = and using 

A = S + &V\ 5$ = 4nSp, and Avi = Here 5 refers to 
an Eulerian perturbation, A refers to a Lagrangian per- 
turbation, and is the Lagrangian displacement. The 
result is the standard perturbed fluid equations with the 
forcing term a* ot (chapter 6 of [6(j): 

^5"-^+^+^=^'. ( 3 - 2 °) 

ot p 



and 



^ = - V 6 . 

p 

Equation (|3.20() can be reexpressed as 



(3.21) 



(3.22) 



where £ is a differential operator [see Eq. (|F2(l ], This 
equation can be solved by expanding £(t,x) in terms of 
a chosen basis of modes £ a (x) and their time-dependent 
amplitudes q a (t), 



€(«,*)= £ g a (t)&»(a:) 



(3.23) 



where a = (n, I, m) label the modes. The basis functions 
can be further expanded in terms of vector spherical har- 
monics (Appendix 0, 



£ a (x) = E£(r)Y E > lm + B£(r)Y B > lm + R<£(r)Y R > lm . 

(3.24) 

The index n — • • • oo is the number of radial nodes, 
and I — • • • oo and m = —I ■ ■ ■ I are the familiar angular 
indices in a spherical harmonic decomposition. We also 
define the inner product and mode normalization 



(€«,£/») = J P(x)& ■ £/3 d 3 X = MR 2 S af: 



(3.25) 



-cj 2 M(x) = jC[£ a (x) 



(3.26) 



Inserting Eq. (|3~23|) into (|X2^1) and using Eq. (|3~2l)|) yields 
the equation of motion for the mode amplitudes, 



MR 2 



(3.27) 



To compute the change in central density we need only 
consider the evolution of the fundamental radial mode 
qo(t)£o(x). [Here and below the subscript refers to 
a = (0,0, 0).] T o see this, substitute Eqs. l{3~2"3"j) and 
(I3~2l|) into iprSIj) . This yields 



Ap 
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^ = £?«(*) 



, Ef R? 

r r 



dRI 



dr 



Y 



(3.28) 

In the r —* limit, Ap/ p must be independent of direc- 
tion (9,6), so it can only be affected by I = rn = radial 
modes 15] . Further, the radial eigenfunctions near the 
center of the star have the form i?^ 00 (r) oc r n+1 , so only 
the fundamental (n — 0) radial mode can change the cen- 
tral density. This radial mode function can be expressed 
as 



(3.29) 



&(aO = if(r)rn = C/(y/4n)$n , 



where C is a normalization constant determined by 
Eq. 13.25f) . n is a unit radial vector, and £g = r [l + 0(r)] 
near r = 0. From Eq. 13.28fl the change in central density 
is 



Sp c 

Pc 



lim — — 

r^O p 



3C 



qo(t) 



(3.30) 



(the Eulerian and Lagrangian density perturbations at 
the center of the star are identical). 

We therefore need to solve Eq. (|3~27|) with a = (0, 0, 0) 
for qo(t). This involves computing the inner product 



(£o,<0 = pr 



C - 



n ■ a tot dft) dr . 



(3.31) 

Using Eqs. lj3~T£)l and (JSHSjl, and the STF integrals in 
Appendix [Bj the angular integral becomes 



n • a tot dQ 



327T 

'"45" 



(3.32) 



The negative sign shows that the angle-averaged radial 
force is inward pointing, leading to compression. 

For a T = 2 polytrope, lu 2 = Att 2 M/R 3 , with 
A « 0.3804, and C w 4.756. The radial integrals 
are computed numerically using the eigenfunction = 
(R/n)£(u) from Appendix |Fj [Eq. l|F10|l ]. giving 



-0.05996M R A B i:j B ij 



(3.33) 
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D. Change in central density 

We now have all the tools needed to compute the 
change in central density at 0(a 7 ). We will specialize 
to a circular binary, for which the tidal fields have the 
form (see Appendix 



£ii £ 



and 



3 sin ujorbt — 2 
-3 cosw or bi sinw or bi 




-3 cos uj OI ht sin uj or ht 

3 COS 2 LUorbt — 2 



Bu = B 



cosworb^ sinw or bt 







COS W or bi 

sincj or bt 




1 

(3.34) 



(3.35) 



KM 



where £ = M'/d 3 , B = 3(M' /d 3 )V olh , w orb 
M')/^ 3 ] 1 / 2 is the Keplerian orbit angular velocity, and 
Vorb = Woi-bd is the relative orbital velocity. Note that 
the tidal fields contracted with themselves do not depend 
on the orbital phase: £ l:j £ lj = 6(M'/d 3 ) 2 and B lJ B l] = 
18(M'/d 3 ) 2 (M+M') jd. Since d evolves very slowly com- 
pared to the orbital and stellar oscillation frequencies, we 
can ignore its time dependence when solving Eq. (|3.27|) 
for qo(t). The initial conditions that qo, qo, and Bij vanish 
at very early times (when the binary is widely separated) 
yield the s imple solut ion q (t) = (£ , a tot ) /(co 2 MR 2 ). 
Equations (|3.3U|I and l|3.33|l then yield the gravitomag- 
netic contribution to the change in central density 



fin R 5 

= 0.06427— BnEV 



(3.36) 



This equation is valid for any slowly-varying magnetic- 
type tidal field, not just the specific form given above. 8 
The formula (|3.35() for Bij along with Eq. (|3.5|l then gives 
the total change in central density up to order 0(a ), 



Sp c 

Pc 



= -2.280 



M'\ 2 




M ) 


(f 




M' 


M 



M' 
~M 



2 T ? 



This formula shows that there is a critical orbital sep- 
aration, d C rit/i? = 0.5074(M/i?) 2 (l + M'/M), where 
the gravitomagnetic crushing force can overwhelm the 
Newtonian tidal stabilization. This separation can be 
large if one considers not only NS/NS binaries but also 
NS/massive BH binaries. However, one must compare 
this separation with an estimate for the onset of tidal 
disruption, or, in the case of massive BHs, the separa- 
tion when the inner-most stable circular orbit (ISCO) 



As another example, consider the magnetic-type tidal field of a 
spinning black hole with spin parameter a = a/M' . Far fro m th e 
hole Bij Bij = 18a 2 M' 4 /d s [from Eqs. (3.36) and (5.45b) of (ril l, 
and the resulting change in central density is still compressive, 
with magnitude Sp c /p c = 1.157a 2 (M / R) 3 (M' /M) 4 (R/d) s . 




FIG. 5: (color online) Critical orbital separation for compres- 
sion to overwhelm stabilization. The solid (blue) curve shows, 
as a function of mass ratio, the critical orbital separation dcrit 
(in units of the NS radius) where the compression and sta- 
bilization terms in Eq. 13. 3711 are equal. The short-dashed 
(green) curve shows the tidal disruption limit dtidai/^?- The 
dashed-dotted (red) curve shows the ISCO di sco /R while the 
long-dashed (black) curve shows the event horizon dhorizon/-R 
of a nonspinning black hole with mass M'. This plot shows 
that in an inspiralling binary, either the tidal disruption limit, 
the ISCO, or the event horizon is reached before the crit- 
ical separation where compression dominates. This is also 
true for rapidly spinning black holes: as the spin parameter 
a/M' varies, the lines corresponding to the ISCO and horizon 
would shift up or down, but they would always remain above 
the curve for d cr it- These curves assume a neutron star with 
M = 1.4Mq, R — 10 km, and slT — 2 polytropic equation of 
state. 



or event horizon is reached (see Fig. [S]). The tidal 
disruption radius is approximately dtidai/-R = 2.4(1 + 
A/'/M)^/ 3 ) HI. For equatorial orbits the ISCO oc- 
curs at a separation of di sco /R = f3i sco (a)(M/R)(M' /M), 
while the event horizon is at a separation of dhorizon/-R = 
Aiorizon (a) (M/R) (M'/M) (these formulas are strictly 
valid only in the limit where M'/M 3> 1, but we apply 
them for all mass ratios). Here /3i SCO (a) and /?horizon(a) 
are dimcnsionless functions of the BH spin parameter a 
and vary from 1 to 9 and 1 to 6, respectively [fj^]. When 
one compares these critical separations to the onset of 
crushing, one finds that, for any plausible value of com- 
paction (M/R < 1) or mass ratio (M'/M > 1), either 
the tidal disruption, ISCO, or horizon radius is reached 
before the gravitomagnetic crushing force dominates over 
tidal stabilization. Therefore, an increase in central den- 
sity by this mechanism cannot occur. 

The effects of gravitomagnetic compression on the sta- 
bility of relativistic NSs could be treated by applying 
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the second-order perturbation methods of this section to 
an initially static and spherically symmetric relativistic 
star. A simpler approach would be to modify Thome's 
[l7| "local-asymptotic-rest-frame" analysis. This modifi- 
cation would supplement Thome's potential energy func- 
tion for the star [his Eq. (7)] with the following terms: 
(1) the gravitomagnetic contribution to the rate of tidal 
work performed on the NS by a distant tidal field, 9 

dW 1 dlu 2„ Si, 

— = _ -Bu-^- ; 3.38 

dt 2 13 dt 3 J dt ' v ; 

and (2) current-quadrupole contributions to the star's 
internal energy. Extending Thome's analysis to gravit- 
omagnetic interactions confirms the main result of this 
section: the weaker gravitomagnetic compression cannot 
overwhelm the electric-type tidal stabilization of an ini- 
tially static neutron star. 

A possible exception is a situation in which a magnetic- 
type tidal field is present but the electric-type tidal field 
is absent or very small. A (somewhat contrived) exam- 
ple of such a situation would be a judicious arrangement 
of at least two less massive bodies (planets) orbiting a 
star, with their masses, distances, inclinations, and or- 
bital phases carefully chosen, leading to a nearly van- 
ishing £ij but nonzero Bij. Although it is very unlikely 
that such situations exist in nature, they show that stars 
can in principle undergo a net compression due to tidal 
effects. 



IV. CHANGE IN CENTRAL DENSITY FROM A 
PREEXISTING VELOCITY FIELD 

In the previous section we considered the change in 
central density caused by a current-quadrupole moment 
that is induced by an external tidal field. Now we con- 
sider the case in which a current-quadrupole moment or 
other velocity field is preexisting in the star rather than 
induced by tidal interactions. This velocity field then 
couples to the external tidal field to change the central 
density. As discussed in Sec. viscosity will damp most 
astrophysical velocity perturbations, but a velocity field 
might arise as an artifact of a numerical simulation. 10 
In the WMM simulations 0, E3, 13 it is possible that 



9 Equation 13.381 was first derived by Zhang lO'fl . The electric- type 
contribution to the tidal work [the first term on the right-hand- 
side of 13.381 1 was s how n to be gau ge and energy-localization 
inva riant by Purdue l64l . Favata l65l . and Booth and Creighton 
lOfl . Although it has not been explicitly calculated, these prop- 
erties should also hold for the gravitomagnetic term in 13.381 . 

r- m odes driv en unstabl e by radiation reaction in rotating hot 
neutron stars 67, 68, 69] are an additional source of a preexisting 
current quadrupole, but their magnitudes are too small to be 
of interest here: For an r-mode in a star with angular speed 
f2, the characteristic velocity of the current quadrupole Sv ~ 
Sij/(MR 2 ) is approximately Sv ~ &RU ~ &(M/ R) 1 / 2 (n/tt c ), 
where Q c = (M/i? 3 ) 1 / 2 and a C 1 parameterizes the r-mode 
amplitude. Also, the rotation of the star provides additional 
support against collapse. 



a current-quadrupole moment arises in the formulation 
of the initial data. In those simulations, two initially- 
corotating single NS solutions are placed on the com- 
putational grid and are allowed to relax (using artificial 
viscosity) to a two-body equilibrium state which is nei- 
ther corotating nor irrotational. Indications of a current 
quadrupolar velocity pattern can be seen in the original 
WMM simulations (see Figure 4b of [T(| and associated 
discussion). However it is not clear how that current 
quadrupole is generated. 

Begin by considering an approximately spherical, non- 
rotating star that satisfies the fluid equations l|3.1|l aug- 
mented by the 1PN external acceleration l|3.2|l . and that 
contains a preexisting velocity Vq . This velocity field can 
generally be expressed as a sum over vector harmonics as 
in Eq. (l3~m 

vo(t,x) = ^2E l v m (t,r)Y E ' lm + B l ™(t,r)Y B ' lm 

I in 

+R l ™(t,r)Y R < lm . (4.1) 

In isolation, Eqs. (|3.1|) are satisfied with a° xt = and 
v = Vq. If we further impose the condition that the 
density is time-independent in the absence of tidal fields, 
the velocity field must satisfy V • (pv ) =0. If we assume 
that the background density p is spherically symmetric, 
p = p(r), then Vq must also satisfy n ■ Vq = V ■ Vq = 
0; it is therefore proportional to a magnetic- type tidal 
field, v = Y,irn Bl v mYB lm - The magnitude of v is not 
known, but we will assume that it is small enough to 
satisfy Vq <C M/R. Since the (vq ■ V)t>o term is small 
in this approximation, Vq ss 0, and the structure of the 
star in isolation is adequately described by the ordinary 
equations of hydrostatic equilibrium [Eq. (|3.6|) ]. This also 
implies that B 1 ™^, r) is independent of t. Assuming that 
order 0(vq) terms are small allows us to neglect various 
1PN terms in the hydrodynamics equations. Other 1PN 
terms are dropped for the reasons discussed in Appendix 

El 

Now allow the external tidal fields to perturb this star. 
Since the background velocity Vq negligibly affects the 
structure of the star in isolation, Lagrangian perturba- 
tions of the star are described by Eqs. i|3.20[l and i|3.21[l . 
The methods used in Sec. 1111 Cl ean then be used to com- 
pute the change in central density. The main step is 
to compute the fundamental radial mode evolution via 
Eq. |3~27|l [with a = (0,0,0)], using (£ ,a cxt ) for the in- 
ner product. Since Vq is small, we ignore terms of order 
0{vl) in af xt . 

Substituting the expansion (|4.1|) for vq into af xt , ex- 
pressing the vector spherical harmonics in terms of STF-Z 
tensors [Eqs. 1B8(| ]. and performing the angular integra- 
tion, we find that the only piece of the velocity field that 
can change the central density is proportional to an I = 2 
magnetic-type vector harmonic, Y B ' 2m , which couples to 
the VqxB piece of the external acceleration. 11 The result 



11 There is also a contribution to the inner product (fofs),^ 1 ') 
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FIG. 6: (color online) Change in central density of a neu- 
tron star with a preexisting velocity field [Eq. (14.611 with the 
cosine term set to +1]. The size of the current quadrupole 
is fixed at Vs = 0.1(A//i?) 1/2 . All curves are for a NS 
with mass M = IAMq and radius R = 10 km. From 
left to right, the first two curves are for a NS compan- 
ion with M'/M — l(blue,solid), 3(red,short-dashed) and are 
terminated at the tidal disruption radius. The next two 
(black) curves are for a black hole companion with M'/M = 
5 (dot-dashed), lO(long-dashed) and are terminated at the 
inner-most stable circular orbit. For the value of Vs used 
here, compression dominates over stabilization before tidal 
disruption or orbit instability occurs. 



FIG. 7: (color online) Same as Figure except that the 
size of the current-quadrupole moment is smaller, Vs = 
0.01(M/R) 1/2 . The compression is much smaller in this case 
and, for three of the stars, is eventually dominated by the 
tidal stabilization. For the M'/M — 3 binary the neutron 
star central density decreases by < 0.5% before tidal disrup- 
tion. 



If we approximate £o( r ) = r [l + 0(r)] w r in Eq. I|4.2|) 
(incurring an error < 20% for a T = 2 polytrope) and 
combine with l|4.3ll . the inner product simplifies to 



for the inner product is 



C 



(4.4) 



C Air 2 



E 



yfr pr*e (r)Bi m (t,r)dr 



(4.2) 



A current-quadrupole moment is also proportional to 
Y B,2m _ substituting Eq. (gj) into (JT3J yields 

Sij= -T\ll £ (Jpr'B^r^yi™. (4.3) 



from a velocity component proportional to y K,2m coupling with 
the 3vi<£ oxt piece of the external acceleration; but this piece is 
excluded by our condition that V ■ (pvg) = when the star is 
in isolation. If we expand the gravitational potentials to higher 
powers of a (including octupole and higher tidal moments), other 
velocity couplings that could change the central density are pos- 
sible, but would be smaller in magnitude. 



Since i)o ~ we can assume that Sij is a constant 
and parameterize the magnitude of its components by 
\Sij\ ~ MR 2 Vs, where V$ is the characteristic velocity 
associated with the current quadrupolar motions. Inte- 
grating Eq. (|3.27l) using (|4.4H and 13. 3511 . assuming that 
d varies slowly compared with the orbital and stellar os- 
cillation periods, and using the condition that <7o> 9o ~~ * 
as t — > — oo, we get 



<Zo (t) 



6C M' /M + M'\ 1/2 V s cos(w orb t + 5) 



where luq is the angular frequency of the fundamen- 
tal radial mode, V s = {S% z + S^f/ 2 /{MR 2 ), and 
tan 5 = —Sy Z / S xz . Since the fundamental mode fre- 
quency ujo/2it 4.2 kHz (for a canonical M — IAMq, 
R = 10 km NS with a T = 2 EOS) is several times larger 
than the orbital frequency (ss 590 Hz at d/R ~ 3 for two 



NSs), we can usually approximate Uq - 



Using 



Eqs. (|3.30|) and 13.5(1 . the total change in central density 



14 



is 



dp 



M'\ 2 /i? x 6 



17.2(i\5 ( — 



1/2 



Because of the cosine term in (14 .6f) , the change in central 
density oscillates in sign. Compression or stabilization 
depends on the orbital phase. Since we are interested in 
the possibility of crushing forces and how they compare 
with the Newtonian tidal stabilization, we will set the 
cosine term to +1 in our discussion below and in Figures 
El and □ 

In contrast to the case treated in Sec. IIIII when the 
current-quadrupole moment is preexisting the gravito- 
magnetic crushing contributes to the change in central 
density with a lower power of a than the stabilizing term. 
This means that at large values of d/R — I /a, crushing 
will dominate over stabilization even if Vs is small. In 
Fig. we plot the total change in central density for a 
1.4M , 10 km NS with V s = 0.1{M/R)^ 2 in a binary 
with mass ratios of M' /M = 1, 3, 5, and 10. The gravito- 
magnetic term clearly dominates, leading to compression. 
If the size of the current quadrupole is reduced by a factor 
of 10 to Vs = 0.01(M/i?) 1/2 , the change in central den- 
sity becomes much smaller (Fig. [7J . While compression 
still dominates at large separations, the tidal stabilization 
eventually overwhelms the gravitomagnetic compression. 
For the stars treated here significant compression would 
require rather large current-quadrupole moments, with 
V s > O.l(MfR) 1 / 2 . 

Although the changes in central density shown in Fig- 
ures and \7\ are small, the gravitomagnetic crushing 
could be enhanced by an orbital resonance with the fun- 
damental mode. Although we have made the approxi- 
mation that ujq ^> ^ rb , this is generally true only for 
Newtonian stars, which do not have a maximum mass. 
For relativistic stars, the fundamental mode frequency 
approaches zero as the mass of the star approaches the 
maximum mass of its EOS. This means that for stars 
sufficiently close to their maximum mass, the fundamen- 
tal frequency could be low enough to be in resonance 
with the orbital period before tidal disruption occurs. 
This resonance would amplify the change in central den- 
sity. Even if a resonance does not occur, a relativistic 
star that is close to its maximum mass is more easily 
perturbed past the critical point of its potential for ra- 
dial oscillations (see Thorne 0]). ^ compressed enough 
such stars could undergo gravitational collapse to BHs — 
although their masses would have to be very close to the 
maximum mass for this to happen. 



V. DISCUSSION AND CONCLUSIONS 

Is the binary-induced compression and collapse of a 
neutron star possible? In the simplest and most realis- 
tic case of two nonrotating, initially unperturbed stars 



in a binary, the answer is no. Although there is a com- 
pressional force on the star it is always smaller than the 
Newtonian force that stabilizes the star. The compres- 
sional force is small because it arises through a nonlinear 
fluid interaction with a post-Newtonian tidal field: the 
gravitomagnetic field induces a current quadrupolar ve- 
locity field which then couples to itself and to the grav- 
itomagnetic field to produce compression. The Newto- 
nian stabilization force also arises through a nonlinear 
fluid interaction, but it is induced by a Newtonian tidal 
field instead. This stabilization force is also small, but 
not as small as the post-Newtonian compressional force. 
Although we only consider first post-Newtonian effects, 
it seems unlikely that effects at 2PN and higher orders 
will be large enough to change the sign of the change 
in central density. For a NS/NS binary the parameters 
e = M/R ~ 0.2 and a = R/d < 0.3 are sufficiently small 
that higher order terms in an expansion of the central 
density in e and a are unlikely to be important (unless the 
coefficients of those terms are much larger than unity). 

However, there are certain physical situations in which 
a net compression is possible in principle. For configura- 
tions of masses in which the Newtonian tidal field nearly 
cancels, the gravitomagnetic compression dominates over 
tidal stabilization. Such configurations probably do not 
occur very frequently in nature. A somewhat more likely 
possibility arises when a velocity field is already present 
in the neutron star and does not need to be induced 
by tidal interactions. If the velocity field has a current 
quadrupolar component, a net compression due to grav- 
itomagnetic forces is possible. In Newtonian stars this 
compression is small for plausible values of the internal 
fluid velocity. Nevertheless, stars that are close to their 
maximum mass could, in principle, be pushed beyond 
their stability limit and made to collapse to black holes. 

The implications of our results for the revised Wilson- 
Mathews ^ij simulations are unclear. As discussed in 
Sec. Ill Al other groups appear to have ruled out com- 
pression in NS /NS simulations that enforce irrotation. It 
therefore seems very likely that the compression seen in 
the irrotational simulations of Marronetti, Mathews, and 
Wilson 01 is unphysical and possibly related to their 
method of implementing certain boundary conditions, or 
is due to insufficient resolution. If low resolution is the 
cause of the compression in their irrotational simulations, 
then it may also be the source of the small residual com- 
pression in their unconstrained hydrodynamics simula- 
tions [l4|. If low resolution is not the source (as main- 
tained by Wilson |z3), then there remains the possibility 
that the compression is caused by Wilson and Mathews' 
method of determining the initial data or the use of arti- 
ficial viscosity. Our analysis in Sec. IIVI shows that com- 
pression can arise if the initial velocity configuration has a 
current quadrupolar component. For current quadrupo- 
lar velocities of size Vs ~ 0.1(M / 'R) 1 / 2 we predict central 
compressions of order Sp c /p c < 1%. This is roughly a fac- 
tor of ~ 5 — 80 times smaller than the compression seen 
in Wilson and Mathews' revised simulations [l4|. Us- 
ing a larger value of Vs could bring our estimates closer 
to the Wilson-Mathews' values, but would begin to vio- 
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late our assumption that Vs is small. The actual size of 
the current-quadrupole moment in the revised Wilson- 
Mathews simulations is not clear. If it is very small, 
then our approach might not be a plausible explanation 
for compression. Higher values for the central compres- 
sion (possibly leading to gravitational collapse) could be 
attained by extending our calculations to include rela- 
tivistic stars near their maximum mass or the effects of 
orbital resonances. While the scenario discussed here is 
plausible in the context of binary neutron star simula- 
tions, in nature a preexisting current quadrupole will be 
rapidly damped and so will be irrelevant for observations 
unless it is generated shortly before coalescence. 

The claims made in this paper should be testable in a 
numerical simulation. Using a binary neutron star simu- 
lation in full general relativity, one could artificially im- 
pose a current quadrupolar velocity field on the stars 
and see if compression results. In a simulation with no 
initial current-quadrupole moment, it should be possible 
to extract information about the gravitomagnetic Love 
number by decomposing the late-time velocity field into 
vector spherical harmonics. Our claims could also be 
verified using simpler Newtonian and relativistic codes 
that model the hydrodynamics of a single neutron star. 
Such codes have been useful in studying the nonlinear 
evolution of r- modes 0, [T^, El] ■ These codes could pre- 
sumably be modified to treat binary systems by adding 
tidal acceleration terms fas in Eq. JsOJ] to the hydrody- 
namics equations. This would also allow the effects of 
magnetic-type and electric-type tidal interactions to be 
tested separately by turning those specific terms on or 
off, something that is harder to do in fully relativistic 
binary NS simulations. 

The revised Wilson-Mathews simulations 01 were per- 
formed over five years ago. It would be interesting to re- 
examine those simulations using higher resolutions than 
were possible at that time. If low resolution or some other 
computational artifact is not the cause of the Wilson- 
Mathews compression, definitively determining the com- 
pression's source will require isolating those pieces of 
physics that are contained in the Wilson-Mathews simu- 
lations but are not contained in other NS/NS codes (none 
of which indicate compression). The Wilson-Mathews 
method of choosing the initial data (which is neither coro- 
tational nor irrotational) and their soft equation of state 
appear to be two areas that are worth further investi- 
gation. In particular further work on generalizing the 
range of initial data sets for NS/NS binaries is encour- 
aged. Real neutron stars are likely to have some spin and 
could be differentially rotating. Such stars could not be 
modelled in simulations that constrain the initial data to 
be corotating or irrotational. Restricting the initial data 
to these two classes could neglect potentially observable 
effects such as spin-interactions and r-mode excitation 

As this paper was nearing completion, a new analysis 
of the compression effect by Miller [74| appeared. Miller 
performed binary NS simulations in full general relativ- 
ity and found that the central density decreases with de- 
creasing orbital separation according to Sp c /p c oc a 1A . 



Because the NSs in his simulations were initially coro- 
tating, Miller's results do not directly contradict those 
of Wilson and Mathews (their stars relax to a configu- 
ration of almost no intrinsic spin; see Sec. [Hj. Further- 
more, Miller's discrepancy with the Sp c /p c cx a 6 scal- 
ing for tidal stabilization is not necessarily surprising 
since calculations of the a 6 scaling 0, 0, 0] assume 
that the neutron stars are nonrotating. It is interest- 
ing to compare Miller's results with the central density 
scalings shown in Fig. 2 of Taniguchi and Gourgoulhon: 
A rough fit to those curves shows that 6p c /p c oc a 5 - 7 
for their irrotational simulation (with values of order 
\Sp c /p c \ ~ 0.002 - 0.01) and 5p c /p c oc a 3 for their 
corotating simulation (with values of order \5p c /p c \ ~ 
0.02 — 0.1). This indicates that the smaller central den- 
sity changes in irrotational simulations are dominated by 
tidal-stabilization effects scaling as Sp c /p c oc a 6 , while 
the larger central density changes in corotating simu- 
lations are dominated by rotational stabilization effects 
scaling as Sp c /p c oc fl 2 oc a 3 [where fl is the orbital 
and rotational angular frequency; see Eq. l|G6|l ]. Al- 
though probably coincidental, it is also interesting to note 
that the revised Wilson-Mathews simulations also found 
a 6p c /p c oc a 1 ' 4 scaling, although with the sign appropri- 
ate for compression (see Table IV of [14] or our Figure 
Miller's Sp c /p c oc a scaling remains puzzling, but 
we speculate that it may arise from his use of initially- 
corotating neutron stars. Examining the central density 
scaling in fully-relativistic simulations of initially irrota- 
tional or arbitrarily spinning neutron stars could help to 
resolve the issues discussed here. 
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APPENDIX A: JUSTIFICATION OF METRIC 
EXPANSION AND EQUATIONS OF MOTION 

The purpose of this appendix is to justify the form of 
the metric in Eq. (|l.l|l and the neglect of certain 1PN 
terms in the metric and hydrodynamics equations (|3.1|) . 
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1. Metric expansion 

The purpose of our metric expansion (|1.1|) is to pro- 
vide a coordinate system in which the properties of a 
star (specifically, the change in central density) interact- 
ing with a binary companion can be studied. In this co- 
ordinate system the binary companion (labelled B here 
and denoted with a prime in the main text) interacts with 
the star (labelled A here) only through tidal interactions. 
Since we wish to study only the leading-order corrections 
to the change in central density due to post-Newtonian 
interactions, we throw away certain 1PN terms; this is 
justified below. Our approximation amounts to keeping 
only the 1PN gravitomagnetic tidal corrections to the 
metric. The justification for the form of our metric will 
rely heavi ly o n the formalism introduced in Racine and 
Flanagan [75|. The tidal pieces of our metric are merely 
the I = 2 tidal pieces of the Newtonian and gravitomag- 
netic parts of the "body-adapted frame" metric derived 
there. We specialize the general treatment given in [7{| 
to the limited context of a Newtonian star interacting 
with quadrupolar tidal fields. The reader is referred to 
that paper for further details. 

Begin by considering the 1PN expansion of the metric 
in terms of the standard PN parameter i = 1/c: 



ds 2 = - {l + 2e 2 <f> g {eT,X) + 2e 4 [<S> 2 g (eT 1 X) 
+ip g {iT, X)} + 0(s 6 )}dT 2 



(Al) 



+ [2e d Cf (eT, X) + 0{e b )]dTdX l 

+ [S i:j - 2e 2 <S> g (sT, X)5 l3 + 0{e i )]dX l dX j . 

This metric describes the global coordinate system of the 
binary. The coordinate system is conformally Cartesian 
and asymptotically flat (since we specify that the global 
potentials $ g , , and ip g go to zero far from the binary). 
We also assume global harmonic gauge, 



da 



d(iT) dX i 



(A2) 



(Note that our notation differs from that used in [73 and 
that our time coordinates differ from theirs by a factor of 
1 /i. Time derivatives of a quantity introduce additional 
factors of e. See Sec. II of |75|.1 

In the vicinity of star A, there exist local coordinate 
systems (t, x) in which the metric has an expansion of 
the same form as (Al), 



ds 2 



{1 + 2e 2 $ loc (et, x) + 2e 4 [$f oc (et, x) 
+^ oc (et,x)} + 0(e 6 )}dt 2 
[2e 3 ( loc (et,x) + 0(e 5 )}dtd. 



(A3) 



x 



and ^ioc diverge at large distances from the star due to 
the tidal contribution to the potentials. 

The local and global coordinate frames are related by 
a 1PN coordinate transformation of the form 

T(t, x j ) = t + ia(it, x r ) + £ 3 (3{it, x j ) + 0{e 5 ) , (A4a) 



X l (t, x 3 ) = x l + z l (et) + e 2 ^(et, x j ) + 0(e 4 ) , (A4b) 

where z l is the Newtonian order spatial vector that re- 
lates the global frame to the local frame. 12 The standard 
coordinate transformation of the metric components, 
along with the gauge conditions (|A2|) in both frames, 
relate the potentials in the global and local frames and 
provide the functional form of a, (3, and h l up to several 
freely specifiable functions of time and one solution of 
Laplace's equation (see Sec. IIB of [zD])- In the vacuum 
region outside the star (but far from the binary com- 
panion), the local potentials ^loc, Cl° c > an d ^loc can be 
expressed as a multipole expansion in powers of r l and 
l/r l+1 , where r — (xiX 1 ) 1 ^ 2 and I is the angu lar harmonic 
index of the expansion [see Eq. (3.28) of |Z3]. These mul- 
tipole expansions are characterized in terms of body mo- 
ments (the coefficients in front of the l/r' +1 terms), tidal 
moments (the coefficients in front of the r l terms), and 
gauge moments (coefficients that appear in front of both 
types of terms and contain information about the coor- 
dinate system, but which do not contain gauge-invariant 
information about the stars). Racine and Flanagan [75| 
show that the freely-specifiable pieces of the functions 
that appear in the coordinate transformation (|A4|) can 
be chosen in such a way that (i) all the gauge moments 
vanish, (ii) the full 1PN mass dipole moment of the star 
vanishes, and (iii) all the tidal moments with I < 2 van- 
ish (see Table I of [75| and the associated discussion). 
These choices uniquely specify a body-adapted harmonic 
coordinate system. 

We further restrict the metric l|A3(l by throwing away 
the nonlinear 1PN terms £ 4 (& 2 oc + if>\oc)- We argue below 
that these terms will not affect the central density at the 
order in which we are interested. Since the remaining po- 
tentials satisfy linear partial differential equations, they 
can be unambiguously split into terms that depend on 
the self-gravity of the star and terms that depend on the 
external tidal field, 



<I> 



loc 



self 



<j> c 



*loc 



cr lf + ci 



+ [Sij - 2e 2 $i oc (£t, x)5 lJ + 0(e 4 )}dx l dx 3 



(A5) 



(A6) 



and in which the gauge condition also has the same form 
as in lA2p . The potentials in both coordinate systems 
satisfy the standard, 1PN Einstein equations in harmonic 
gauge. However, the metric in the local frame of the star 
is not asymptotically flat as the local potentials ^loc, C]° C ! 



12 Terms in the coordinate transformation at higher order in i do 
not affect the metric at fPN order. Also, terms at 0(e 2 ) in 
T(t,xi) and 0(i) in Jf'Jt,! 3 ) produce only constant shifts of 
the coordinate systems and can be set to zero. 
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(In the body of this paper 

$seif = $ ) The externa i 
potentials satisfy the vacuum Einstein equations 

V 2 $ cxt = V 2£cxt = q ( A7 ) 

and are given by the r l pieces of the multipole expansions 
of $ loc and C' oc - In Racine and Flanagan jf^, the self 
pieces of the potentials also satisfy the vacuum equations 
as in (| AT|) and are given by the l/r l+1 pieces of the multi- 
pole expansion. Such an expansion diverges at the center 
of the coordinate system (when r — > 0). Since they are 
interested in treating the dynamics of strongly gravitat- 
ing bodies, the "body-adapted" coordinates of are 
only valid in the weak-field "buffer region" that exists 
outside the body but far from the companion. In this 
paper we wish to treat the internal dynamics of a weakly 
gravitating fluid star. To do this we use the slightly modi- 
fied body-adapted frame described in the last paragraph 
of Sec. Ill D of [73, which extends smoothly into the 
star's interior. The self potentials are given by the usual 
Poisson integrals associated with the equations 

V 2 $ sclf = 47rp and (A8) 



V 2 Cf lf = , (A9) 

where p is the mass density and Vi is the fluid velocity. 

Since we are interested in the oscillation modes of the 
star, the explicit form of the metric outside the star does 
not concern us. To find the explicit form of the self po- 
tentials inside star, one would simply solve Equations 
(|A8|) and (|A9|> along with the hydrodynamic equations. 
However, for our purposes we can also ignore the Q 
piece of the metric as it will not affect our calculation of 
the change in central density; this is justified below. To 
explicitly compute the external potentials (both inside 
and outside the star), one can use the formulas found in 
Sec. V of [75j. Computing only the lowest order (I = 2) 
piece of those potentials yields 

$ oxt = -Eijx'a? (A10) 

and 

C- Xt = -\Y ijk x i x k . (All) 

In the notation of |2a|, £ij = — n Gfj and, for a two-body 
system, is given by 

M B ( zf A zf A \ 
E l3 = (% - ^ , (A12) 

where Mb is the mass of the companion {M 1 in the body 
of this paper) , and zf A is the separation vector pointing 
from body B to body A. The tidal moment Yy/- is given 
by 

Yijk — Aijk — A<ijk> , (A13) 



where 

A ijk = -AS jk Vf A + I {Sij'zi + S ik zf) - U jk z A , 

(A14) 

and <> means to symmetrize and remove traces on the 
enclosed indices. Using the definitions By = — |V \jB^ 
and B = V x C° xt ; the magnetic-type tidal field By can 
be related to the Yij k tidal moments by 

®ij = j( e iabY baJ + ejabYbai) , (A15) 

and the external gravitomagnetic potential can be ex- 
pressed as 

Cr* = -\e iaj B a k x^x k . (A16) 

Racine and Flanagan |58j give an equivalent but simpler 
formula for the tidal field: 

B V} = GMsz^e^z^V^/d 5 . (A17) 

In the above equations, zf is the position vector of 
body A relative to the center of mass of the binary, 
yBA _ 2.BA j g ^ ne r el a tive velocity of the bodies, zf A = 
zf — zf, d = \zf A \, and zf = M B zf A jd z . For a circular, 
Newtonian orbit in the x-y plane, zf A has the Cartesian 
components 

zf A = d[cos cj or bt, sin w or b*, 0] , (A18) 

where w orb = [{M A + M B )/d 3 ] 1 / 2 and zf = 
—MbzP a /(Ma + Mb)- In this case the components of 
£ij and By are given by Eqs. (|3.34|) and (|3.35|l . 

We note that the expression (|A10|I and (|A16|I for the 
external pieces of the metric are identical to standard 
expressions for the metric in the vicinity of a point parti- 
cle in an arbitrary gravitational field expressed in Fermi 
normal coordinates [lilll]}- In this language, consider 
the worldline z a {t) of an observer moving on a geodesic 
near a gravitating body. Manasse and Misner [76| have 
shown that one can introduce a coordinate system in 
which g^ v = rj^ and r M „ CT = are satisfied all along this 
worldline. Such coordinate systems are called Fermi nor- 
mal coordinates and describe the proper reference frame 
of a freely falling observer. Near the origin of this coordi- 
nate system in a specific choice of gauge, the metric can 
be expanded as a power series in the distance la^l from 
the observer as [33, UM 

ds 2 = (— 1 — RoiojX l x j )dt 2 + ^RoijkX^x 3 J dtdx k 

+ (s kl - ^R lk] ix l x^ dx k dx l + 0(\x j \ 3 )dx a dx !3 , (A19) 

where the coordinate time t is also the proper time 
along the observer's worldline, and R a /3-yS are compo- 
nents of the Riemann tensor evaluated along the geodesic 
worldline z a (t). Some of these components can be ex- 
pressed in terms of the tidal moments £y (t) = Roioj 
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and Bij(t) = ^ei qp RQj pq . These moments describe the 
tidal field of a strongly relativistic object, but in the 
PN (weak-field) limit they reduce to the tidal moments 
described above [Eqs. I|A12|) and (|A15|l ]. Higher order 
tidal moments are defined in terms of the derivatives of 
the Riemann tensor evaluated along the worldline. Ishii 
et al. [73 have recently extended the metric expansion 
(|A19|) to Ofla^'l )■ For the extension to accelerated and 
rotatin g ob servers, see Ni and Zimmermann |78( and Li 

and Ni [23. 

The post-Newtonian limit of the metric (|A19|) is equiv- 
alent to keeping only the <E> ext and (f xt terms in the 
metric of Eq. jl.ip . To arrive at the full metric (|1.1|) 
for a Newtonian star in a tidal field, we must include 
the Newtonian potential $ sclf . Since we treat the star's 
self-gravity at Newtonian order and neglect all nonlin- 
ear gravitational terms, the superposition principle al- 
lows this extension to be achieved by simply adding the 
appropriate $ sclf terms to the time-time and space-space 
pieces of the metric (|A19|) . The resulting metric and the 
hydrodynamic equations that are derived using it are of- 
ten used in studies of the tidal disruption of an ordinary 
star or compact object [see jT^} and references therein, 
and also their Eq. (127)]. 



2. Neglecting 1PN terms in initially unperturbed 
stars 

We now explain why certain 1PN terms will not affect 
our calculations and can be dropped from the equations 
of motion. For the case in which the star is initially un- 
perturbed fSec. Illlf) . three facts help us to justify drop- 
ping the irrelevant terms: First, since the change in cen- 
tral density is due only to tidal interactions, Sp c /p c must 
be constructed from specific combinations of tidal mo- 
ments. Since these combinations must have even parity, 
the change in central density must have the form 

^ = aSijSV + c 2 BijB ij 

Pc 

+ c 3 £ ijk £ ijk + Ci B ijk 13 ijk + ■■■ , (A20) 

where the coefficients c\, C2, • • ■ depend on M, R, and the 
equation of state. Terms of the form EijB^ are excluded 
because they are parity odd. Since we are interested in 
only the leading-order correction to the Newtonian tidal- 
stabilization term — the B t jB 1 ^ ~ 0(a 7 ) term in l|A20jl — 
it is clear that any 1PN tidal terms that depend on 
electric-type tidal fields £ ai a 2 ---ai can be excluded. This 
immediately excludes all the terms in af xt that depend 
on <E> ext [see Eq. (|3.2|) ]. It also excludes the tidal pieces 

of + 

Second, we are uninterested in O(e) corrections to each 
of the terms in (|A20I) . For example, several terms in 
the 1PN hydrodynamics equations will effect the change 
in central density by adding corrections to the first two 
terms in lfA"20|) that scale like Sp c /p c - [0(1) + 0(e)]a 6 + 
[0(e 2 )+0(e 3 )]a 7 +0(a 8 ). We drop terms that contribute 



to these O(e) and 0(e 3 ) corrections. 13 

Third, any acceleration terms in the equations of mo- 
tion whose angle-averaged radial piece (n ■ a) vanishes 
cannot contribute at linear order to the change in cen- 
tral density. Such terms are also dropped [except in the 
linearized hydrodynamic equations H3.7fl - I|3.9[l ]. All 1PN 
terms are excluded from the metric (|1.1[) or equations of 
motion 1)3. 1|) and i|3.2fl because of one or more of these 
three reasons. 

For example, let us consider some of the terms appear- 
ing in the 1PN Euler's equation in detail. In the absence 
of tidal forces (for a static star), the 1PN terms mod- 
ify the equations of hydrostatic equilibrium [Eqs. (|3 . Of) ] . 
causing O(e) changes to the background structure of 
the star, including the mass and radius. This will lead 
to 0(e) corrections to the mode functions and to the 
leading-order coefficients in (|A20|) . They can therefore 
be excluded. Now consider perturbations to the 1PN 
terms due to tidal interactions. The nonlinear 1PN terms 
$^ oc + -01OC are dropped because their pieces either mod- 
ify only the background star at O(e), or contain only 
electric-type tidal interactions that either lead to 0(e) 
corrections to the 0(a 6 ) piece of 5p c / p c or have vanish- 

. sc if 

ing (n-a). Consider next the acceleration term dj ~ Q 
which is driven by the velocity uj 1 induced in the star. 

This term scales like Q** ~ P 2 ^orbP (0 M^ ~ e 3 a 5 /R 
and causes a change in density Sp/p ~ e 2 a 5 . However, 
the change in central density from this term vanishes be- 
cause (n-£ sclf ) = (n-i/ 1 )) = [since i>W is purely axial; 
see Sec. 1111 H] , The acceleration term v^ 1 ' x (V x £ self ) 

scales like ~ Rp^v^ ~ e 4 a 7 /R and has a change in 
density that scales like Sp/p ~ e 3 a 7 , providing an O(e) 
correction to the 0(a 7 ) compression term. Terms in the 
metric and equations of motion that depend on £f olf can 
therefore be safely neglected. 

Terms in the 1PN Euler's equation [see Eq. (9.8.15) of 
j59| ] proportional to (P/p)Vi$, ($/p)ViP and $Vi$, 
where $ = <I> SGlf , are dropped because they involve 
only electric-type tidal perturbations. Terms that are 
quadratic in the fluid velocity such as ZJ 2 Vi$, i>i(u fc V/ £ )$, 
p^k[vkVi(P - 2p$ + <Fu 2 )], and p- 1 d/dt(pv 2 v i ) add 
corrections to the central density at higher orders than 
concern us (or vanish completely) . Terms that are linear 
in the fluid velocity like p~ 1 d/dt[v i (P — 2p$)] and 
are also either higher order or vanish identically [since 

Except for terms proportional to 4 ext , none of the 1PN 
terms can contribute to the pieces of the central density 
that we seek, so they are safely dropped from the met- 
ric and equations of motion. The resulting equations of 
motion that we use [the continuity equation, Poisson's 



The change in density caused by an acceleration term ~ a in the 
perturbed 1PN Euler equations has the scaling Sp/p ~ ^/R ~ 
a/(Rw 2 ) ~ (R 2 /M)a, where £ is the displacement of the star 
caused by the acceleration a, and lu 2 ~ M/R 3 is the character- 
istic frequency response of the star. 
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equation, and Euler's equation supplemented by gravit- 
omagnetic terms involving £° xt ; Eqs. (|3.1|) and JUS}] are 
identical to the weak-field, I = 2 tidal-order limit of the 
equations used by Ishii et al. [73 in their Fermi normal 
coordinate description of a Newtonian star interacting 
with a black hole tidal field [see Eq. (127) of [73 and the 
surrounding discussion] . 



APPENDIX B: VECTOR SPHERICAL 
HARMONICS AND STF INTEGRALS 

In this appendix we supply useful formulae regarding 
vector spherical harmonics and STF integrals that we use 
throughout this paper. See Sec. II of Thorne [29] for a 
more detailed discussion. 

We begin with the expansion of the scalar spherical 
harmonics in terms of radial unit vectors, 



Y lm {6,<t>) = ytN Al , 



(Bl) 



and the orthonormality relation 

Y lm Y l ' m '*dn = 6m8„ 



(B4) 



The pure-spin vector harmonics are defined by |29j : 



rE, lm 



x x VY lr ' 



= —n x Y 



B, lm 



VHl + T) 



= n x Y 



E. Ir. 



(B5a) 

(B5b) 
(B5c) 

(B6) 



where J = E, B, or R, and their complex conjugates sat- 
isfy 



They satisfy the orthonormality relation 

yJJm . yj'j'm'* dQ = S JJf S w 8 n 



where Ai = {a x a 2 ■ ■ ■ a/), N Al = n ai n a2 ...n ai and n 3 = 
Xj/r. The coefficients of this expansion, y^?, are sym- 
metric trace- free tensors of rank I (STF-/ tensors). An 
explicit formula for the STF-Z tensors y l j^ is given in 
Eq. (2.12) of Thorne [2^. Their explicit components up 
to I = 2 in a Cartesian coordinate basis [e x , e y , e z ] are: 
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y 



j 
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(B2a) 

(B2b) 
(B2c) 
(B2d) 

(B2e) 
(B2f) 



The scalar spherical harmonics satisfy the eigenvalue 
equation 



2\rlm 



V Z Y 



1(1 + 1) 



Y 



lm 



(B3) 



Y 



J, lm* 



{-l) m Y 



m \r J. I — m 



(B7) 



These vector spherical harmonics can also be expressed 
in terms of the STF-/ tensors: 



E, lm 



[yjX-^-A-i - ^yt N M\ , (B8a) 



Y 



B, lm 



j—j €j pq n p y l q T 2 i ^ i NA l _ 1 



Y R > lm = nj y l £N Al 



(B8b) 
3c) 



The gradient, divergence, and curl of the vector har- 
monics are given by the following formulas: 



VjViY lm , 



(B9a) 



B. lm 



R, lm 



djkY k E ' lm + enkmVjY^ lm , (B9b) 



= - [(5ij - n inj )Y l 



+ y/l(l + l)mY. 



E, In 



(B2g) 


^ . i™ 1 - - 


r 






o, 


(B2h) 


^ . y^> ^ m - 


^ \rlm 
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(B2i) 


V x Y E ' lm = 









Y 



lm 



V x Y B ' lm = 



„(V • Y E ' lm ) - -Y 1 
r 



-(n-v)y 



E, lm 



V x Y 



i?, lm 



VKl+T) Y B,lm 



(B9c) 

(BlOa) 
(BlOb) 
(BlOc) 

(Blla) 

(Bllb) 
(Bile) 
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When performing angular integrals over STF-^ tensors, 
the following integrals over products of unit vectors are 
useful: 



j>N A2l+1 dQ = 0, 



(B12a) 



a 2 u a 3 a 4 



<5a 2 i_ia 2 



+ all distinct permutations] . (B12b) 
More explicitly, the first few nonvanishing integrals are: 

(B13a) 



An 

n a n b dtt = —S ab , 



47T 

n a n b n c n d — — {S ab 5 cd + S ac S bd + S ad S bc ) , (B13b) 
15 

n a n b n c n d n e nf dtt = 

y^r (S ab 5 cd 5 e f + S ab S ce S d f + S ab S c fS de 
+ 5 ac 8 bd d e f + 8 ac 8 be 8 d f + 5 ac S b fS de 
+ S ad S bc S e f + S ad S be Scf + S ad S b f5 ce 

+ <5aefibc8 d f + S ae S bd S c f + 5 ae 5 b f5 cd 

+ S af S bc S de + 5 af S bd S ce + S af 6 be 5 cf ) . (B13c) 



We also frequently use the identity 



(B14) 



APPENDIX C: RELATIVISTIC CIRCULATION 
AND THE IRROTATIONAL APPROXIMATION 

The purpose of this appendix is to explain why flu- 
ids that satisfy the relativistic irrotational condition can 
nonetheless have Newtonian vorticity. We also provide 
an alternative derivation of the gravitomagnetically in- 
duced velocity l|3.11|) . 

Begin bv_ considering the relativistic definition of cir- 
culation 



C = (j) hu a d\ a , 



(CI) 



where h is the specific enthalpy, u a is the fluid four- 
velocity, and the integral is taken around a closed space- 
like curve A. Kelvin's circulation theorem states that 
the relativistic circulation is conserved for a curve A 
that moves with the fluid. If we work in the PN limit 
and assume that our closed curve A is purely spatial 
d\ a = (0,dx l ), we can use Stokes' theorem to write the 
circulation as 

C = (f(h Ul )dx l = /(V x hu), ■ dS l , (C2) 
J\ Js 

where u = Uj and dS l is the normal to the surface 
S whose boundary is A. In the relativistic irrotational 
approximation, hu a — V a $, so C = 0. So while the 



irrotational condition in Newtonian theory is simply 14 
uj = V x v = (where v is the coordinate velocity dx / dt) , 
its relativistic generalization is 



h{S7 x u) + Vh x u = 



(C3) 



This indicates that a relativistic, irrotational fluid can 
have nonvanishing Newtonian circulation so long as the 
fluid velocity is restricted by (|C3|) . 

At 1PN order we can write the circulation C more ex- 
plicitly: Using the metric expansion i|A3(l . it, 
vP = u v- 1 , it c 

polytropes) h = 1 + i 2 (e + P/p), we have 
hui = ev l 



1 + i 2 (v 2 - $i oc ) + (9(e 4 ), and (for 



C!° c + ' 



1 



3$ 



loc 



(C4) 
+ 0(e 5 ), 



where v l is the fluid 3-velocity, p is the mass density, 
peo is the internal energy density, and P is the pres- 
sure. Substituting (|C4|) into l|C2|) gives the circulation 
up to 1PN order. Formally, this 1PN expansion assumes 
the scalings v i ~ e ~ (M/R) 1 / 2 ~ (M/d) 1 / 2 . How- 
ever, when the fluid velocity is generated by tidal in- 
teractions, its scaling is actually much smaller: v l ~ 
(M/i?) 1 / 2 (i?/d) 9 / 2 for electric-type tidal interactions and 
v l ~ (M/i?) 3 / 2 (i?/d) 7 / 2 for magnetic-type tidal interac- 
tions. This allows us to ignore all but the first two terms 
on the right-hand-side of (|C4|) as well as the contribution 
from Q , yielding 



C 



f[Vx(v 
Js 



(C5) 



[see also Eq. (8) of Shapiro |3C|]. By Kelvin's circulation 
theorem, an irrotational fluid satisfies C = for all times, 
and V x v = — V x £ cxt . This has the solution v = 
_£ext _|_ for any scalar function A. The boundary 
condition that v and C° xt — > at early times yields the 
solution given in Eq. i|3.11l) . 



APPENDIX D: DEFINITIONS OF NEWTONIAN 
AND GRAVITOMAGNETIC LOVE NUMBERS 

Here we briefly review the definition of the Newto- 
nian Love number &2 which relates the induced mass 
quadrupole moment to the external Newtonian tidal 
field £ij. We also show how the gravitomagnetic Love 
number 72 is related to the equation of state of a star. 

Consider a fluid star with mass M and mean radius 
i?, interacting with the tidal field of a distant companion 



14 A more precise definition of the relativistic irrotational condi- 
tion is that the rel ativ istic vorticity tensor, w M „ = V„(/m M ) — 
V„(/iUi/), vanishes l48l . Using Euler's equation for a perfect fluid 
in the form ji m V (J (/m„) + V„/i = 0, one can show that u^u^v = 
and £^uii_i V = 0, where £g is the Lie derivative along the fluid 4- 
vclocity. This last equation is the differential version of Kelvin's 
circulation theorem. 
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with mass M' located at position x' = (r' = d, 9' , <f)') in 
a coordinate system centered on the first star's center of 
mass. The distant companion tidally deforms the fluid 
star, giving it a mass quadrupole moment Z^. Near the 
star but outside its surface, the tidal expansion of the 
Newtonian potential up to I = 2 is 

$ = fc 2 — j= — ^P 2 (cos 9) . (Dl) 

r er er 

Here P2 is a Legendre polynomial and 9 is the angle 
between x and x' [see Eqs. (4.121) and (4.155) of 31]]. 
The tidal Love number k 2 is defined to be the dimen- 
sionless coefficient of the 1/r 3 piece of the expansion of 
the potential of a tidally deformed, nonrotating body. To 
show how k 2 is related to the STF moments Zjj and £ ^ , 
expand the Newtonian potential as 

M ZTijwJ-ni 1 9 „ • • . . 
$ = "T~2^^ + 2 r n (D2) 

[this is equivalent to Eq. If) expressed in different no- 
tation] . Expanding the Legendre polynomial as 



P 2 (cos9) 



-in 

Y 

4tt 



J2 Y 2m *(9', 4>')Y 2m (6, <f>) (D3) 



m=-2 
2 



£ ylfri' a n' b y^n inj (D4) 



m=-2 



and equating the corresponding 1/r 3 and r 2 terms in 
Eqs. I|D1|) and l|D2|l . it is easily shown that the mass 
quadrupole moment and Newtonian tidal moment are re- 
lated to the Love number k 2 by 



ly — ~k 2 R £ij 



(D5) 



The value of the Love number k 2 depends on the equa- 
tion of state of the star. To determine this dependence, 
one must solve for the structure of the deformed star. 
Lai, Rasio, & Shapiro |81j have done this for the case of 
tidally deformed, compressible ellipsoids. For ellipsoids 
with a polytropic E9S with index n, the moments of 
inertia are 

Iij = J px,iXj d 3 x — —K n MafSij (no sum) , (E>6) 

where = i?(l + oti) are the axes of the ellipsoid, and, 
for a tidal field along the x-axis ^l| > 



5 M' (R 



a 2 = a 3 = ■ 



(D7) 



In these equations EOS information is contained in the 
coefficients 



5 jf 

3 tit\Q'{(ii)\ 



(D8) 



and q n — K n (l — n/5). Here 9(£) is a solution of the Lane- 
Emden equations, £1 is the first root of this solution, and 
9' = d9/dC 

Expanding (|D6() to linear order in a.;, taking the STF 
piece, and using the form for £y due to a point mass on 
the x-axis gives the relation 



J- — — 
y - 2 



R £i 



(D9) 



The Love number is therefore related to the structure of 
the star by 



(D10) 



For a uniform density star, k 2 = 3/2; for a T — 2 poly- 
trope, k 2 = (10/3) (1 - 6/vr 2 ) 2 fa 0.5124. 

For a nonrotating star immersed in a magnetic-type 
tidal field, we have shown in Sec. IIII Bl that a current- 
quadrupole moment is induced and is related to the tidal 
field by 



(Dll) 



Here 72 is the gravitomagnetic Love number analogous to 
k 2 . For a spherical, Newtonian polytrope, 



72 



2 eo n g 
15 ei\o'(ti)\ 



(D12) 



To show this, perform the change of variables p = 
p c 9 n and r = i?£/£i in Eq. (pTTT^I and use p c = 
Ci/(47r|6»'(^i)|)(A//i? 3 ). For a uniform density star, 72 = 
2/35; for a T = 2 polytrope 72 = 2(tt 4 - 20tt 2 + 
120)/(15tt 4 ) « 0.0274. As is the case with the New- 
tonian Love number, the gravitomagnetic Love number 
becomes smaller as the star becomes more centrally con- 
densed. See Poisson |82| for an extension of the concept of 
gravitomagnetic Love numbers to tidally distorted black 
holes. 



APPENDIX E: TIDAL STABILIZATION FOR A 
NONROTATING NEWTONIAN STAR AT 
ORDER 0(a 6 ) 

To compute the change in central density at order 
0(a 6 ), one could follow a perturbative procedure sim- 
ilar to that used in Sec. IIIII above, expanding in S£ (with 
£b = 0). However, the solutions at order 0(e\) are 
more complicated than those at order O(eg). A sim- 
pler method based on an energy variational principle was 
used by Lai |ll| . In this appendix we provide a con- 
cise derivation of the leading-order contribution to the 
change in central density using the method described by 
Lai but specialized to nonrelativistic stars. The re- 
sulting analytic expression for Sp c / p c is not provided in 
|ll| but ap prox imates the analytic result of Taniguchi & 
Nakamura [jjj which was arrived at by a more difficult 
method. 
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The energy E(p c ) of an isolated, nonrotating star 
with baryon mass M, radius R, and polytropic EOS 
P = Kp 1+1 / n , as a function of its central density p c , 
is 



E(p c ) = E int (p c ) + E glav (p c ) 



(El) 



where 



n M 2 

E^t = / udm = k x KMp x J n = — (E2) 



5-n R 



is the star's internal energy, and 



E e 



[ ™dm = -k 2 M 5 / 3 p 1 J 3 = — 

J r ' 5 — n 



3 M 2 
^n~R~ 



(E3) 
■60]) 



is its gravitational potential energy (see chapter 6 of 60] ). 
In the above equations u = nPj p is the internal energy 
per unit mass, m = m(r) is the enclosed mass as a func- 
tion of radial coordinate r (satisfying dm/dr = Airr 2 p), 
and the various constants are given by 



(E4) 



fc 2 = 



47T0'(fi 



£l 



1/3 
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iefrr 



1/3 



= — w 0.8129 . (E5) 



As in Appendix iDl the function is a solution to the 
Lane-Emden equation, £j is the first root of this solution, 
and 6*' = d6/d^. The numerical values here and below are 
for n = 1. (In this section &2 is not to be confused with 
the Newtonian Love number defined in Appendix [D]) 

The equilibrium central density for a stable star, p c .o, 
lies at the stable minimum of E(p c ) — where dE/dp c = 
and d 2 E/dp 2 > 0. Placing the star in an electric-type 
tidal field modifies the total energy to 



E(pc) = E int {p c ) + E grav (p c ) + Wt(pc) 



(E6) 



where Wt accounts for the interaction energy between 
the tidal field and the star's induced mass quadrupole 
moment, as well as the modification to the star's self- 
gravitational potential energy due to the redistribution 
of mass by the tidal field (the kinetic energy of internal 
fluid oscillations is neglected). For an ellipsoidal star 

MM, 



w t = -x 



M' 2 
~R 



R 



M 



5/3 



- A 'I,^J T- (E7) 



where 



4 " V 5 J 

A' = ^2 f x _ n\ 
4 " V 5 J 



3tt 



4 (tt 2 - 6) 2 w 0.2562 , (E8) 
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47r|0'(a) 



5/3 



5(4tt 



5\l/3 



48tt 4 



-(tt 2 -6) 2 w 0.1713 



(E9) 



and 



5jf^^ = 5 
3 3 



1 - 



0.6535 



(E10) 



In the above equations, we have used M/i? 3 = 
47iyj c |0'(£i)|/£i. The central density in the presence 
of a binary companion can then be determined from 
dE/dp c = dE/dp c + dWt/dp c — 0, the roots of which 
must generally be found numerically. Alternatively, one 
can Taylor expand dE/dp c about the density for an iso- 
lated star, 



dE 
dp c 

where Sp c 
density, 



where 



dE 
dp c 



d 2 E 



S Pc + 0(6p 2 c 



(Ell) 



Pc ~ Pc.o- This yields the change in central 

'i (dw t /dp c y 



Sp c 
Pc 



Pc,0 = 77 



p (d 2 E/dp 2 ) 



i i j-2/3 \ 3n/(3— n) 



3 k! K 



(E12) 



(E13) 



The polytropic constant can be expressed in terms of the 
mass and radius by 



K 



4tt 



cr + >'(6)i r 



V« M l-l/n jR 3/n-l 



1 



(E14) 



Since d 2 E/dp 2 and dW t /dp c are both positive, we see 
that at order 0(a 6 ) stars are stabilized by an external 
tidal field. For n = 1, K = 2R 2 /n and Eq. (|ET2|) reduces 
to 



£~£c-o-(£)'(f 



(E15) 



Taniguchi & Nakamura a l so calculated the leading- 
order change in central density. Instead of an energy 
variational principle for an ellipsoidal star, they solved 
the Newtonian fluid equations perturbatively up to order 
0(a 6 ) for irrotational stars. Their result for the change 
in central density is 



Spc 
Pc 



45 

27T 2 



M' 
~M 



(E16) 



which agrees with Eq. I|E15|I to ~ 10%. 



APPENDIX F: FUNDAMENTAL RADIAL MODE 
OF A T = 2 NEWTONIAN POLYTROPE 

Since changes to a star's central density occur along 
the fundamental radial mode of oscillation, we will need 
to compute the frequency of this mode and its corre- 
sponding eigenfunction. We specialize to a Newtonian 
polytrope with EOS P = Kp r and T = 2. 
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For a Y = 2 polytrope the Lane-Emden equations (see 
e.g., chapter 3 of |6(j) yield the following solutions for the 
internal density, pressure, gravitational potential, and 
mass distribution: 



p{u) 



S1I1M 



P(u) 



2 2 p2 sin 2 M 
—p c n — — 

TT U z 



TT 



smu 



m(u) 



jp c R 3 (sin u — u cos u) 



(Fla) 



(Fib) 



(Flc) 



(Fid) 



where u — irr/R and the central density is p c — 
(tt/4)M/R 3 . A T = 2 polytrope also satisfies M/R = 
2Kp c and R = (ttK/2) 1 / 2 . (In this section u is the renor- 
malized radial coordinate and the symbol £ is reserved for 
the mode function.) 

The oscillations of a Newtonian star are described by 
the eigenvalue equation H3.26[) . where 

I'CX = Vi (I V'V.i"} - (V^) ViP 

+ (Vi£ j ) VjP - p&VjVi* - pVi5® .(F2) 

For radial perturbations £ = £(r)n , this equation sim- 
plifies to 



d_ 

dr 



AdP 
r dr 



£ + W 2 p£ = 0, (F3) 



(see chapter 6 of |6(|). In the above equations Y\ is the 
adiabatic index for the perturbations, which we will take 
to be Ti = r = 2 throughout the star. Mode indices are 
ignored throughout this section. For a Y — 2 polytrope 
we can use Eqs. (IF1|) to rewrite <|F3() as 



it 2 sin u £ uu + 2u 2 cos u £ jU 



4 | 1 — — | CD> // 



- ! : > - jT ) siu '< + Au 3 



£ = 0, (F4) 



where £ jU = d^/du and A = oj 2 / (2irYip c ). 

The boundary conditions Eq. (|F3|) must satisfy are 
£(r = 0) = and ^(r = R) is finite. However, more 
specific boundary conditions are needed in order to solve 
the eigenvalue problem. To determine these conditions 
we analytically explore the solutions to Eq. (|F4|> near 
the origin and the stellar surface. Near u — we 
Taylor expand cosm = 1 — u 2 /2 + u 4 /24 — ••• and 
sinw = u(l — u 2 /Q + m 4 /120 — • • • ) and look for a se- 
ries solution of the form £ = X)^Lo a « u " '• Substituting 
into (|F4ll and solving order by order we arrive at the 
recursion relation 



O-n+2 



5n — 6/3 



6(n + 4)(n + l) 



(F5a) 



0,4 







(F5b) 



where (3 = A — 1 + 8/(3Ti). The approximate solution 
near the origin is therefore £(it) ~ a\u + a^u 3 + ■ ■ ■ . To 
find a solution near the surface we use a similar proce- 
dure: Taylor expand l|F4(l about u = ir and look for a 
power series solution of the form £ = J2m=o bm{n — u) m . 
The resulting leading-order terms are £ (u) ~ bo + bi (ir — 
«) + •••, where 



6i 
6o 



1 



2 

r7 



nA 
~2~ 



(F6) 



The constants a± and 6o are undetermined. 

We now outline a numerical "shooting" method to 
compute the eigenvalue A and the eigenfunction (1) 
Pick an initial guess for A. (2) Since the eigenfunctions 
are only defined up to a normalization constant, choose 
a\ = 1. Integrate Eq. I|F4|I as an initial value problem, 
starting a small distance u = 5 from the center of the 
star with initial conditions £(5) =5 and = 1. (3) 

Integrate to u — 7r and compute the boundary condition 



+ (&iAo)£W = o 



(F7) 



Initially this equation will not be satisfied. (4) Choose 
a new value for A such that, when the equations are 
integrated again, Eq. (|F7|) is more closely satisfied. (5) 
Repeat this procedure until l|F7|) is satisfied to the desired 
precision. 

To determine the fundamental radial mode, we want to 
choose a low enough value for A so that the eigenfunction 
has no nodes. Higher frequency radial modes can be 
found by choosing A such that n nodes appear in the 
eigenfunction (where n is the radial quantum number 
for the mode). Following the above procedure we find 
A = 0.3804 for the eigenvalue. Fitting a seventh-order 
polynomial to the eigenfunction gives 

= u- 0.001119m 2 + 0.03302w 3 - 0.006397m 4 
+0.005408W 5 - 0.001502u 6 + 0.0002397m 7 . (F8) 

This function is related to the mode function in Sec. lIII 01 
by £o(x) = (C / \Z4n)(R/TT)£ i (u)n. The normalization 
constant is found from Eq. (|3.25|) : 



C = 2tt 2 



m£ 2 (m) sin m du 



-1/2 



4.756 



In computing the inner product in Eq. 
make use of the integral 



m 4 £(m) sin m du ~ 76.93 



(F9) 



we also 



(F10) 



APPENDIX G: CHANGE IN CENTRAL 
DENSITY FOR ROTATING STARS 

In this appendix, we briefly discuss how one would ex- 
tend our computations in Sec. IIIII to an initially unper- 
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turbed, slowly-rotating star with uniform angular veloc- 
ity r/tt. We show that gravitomagnetic contributions to 
the change in central density vanish at orders 0(e B rj 1 ) 
and 0(e B rj 2 ), but do not necessarily vanish at orders 
Oielrj 1 ) and 0(e 2 B r] 2 ). 

Along the lines of Eqs. I|3.4[l . we expand the fluid vari- 
ables in two dimensionless parameters, rj and e, that char- 
acterize the spin of the star and the external tidal field 
acting on it: 



Solving the analog of Eq. I|3.27|) . and using Eq. H3.30I1 
and the condition that qo, qo —> at t — > — oo yields the 
change in central density at this order, 



,(0,2) 



(t,0) 



p(°.°)(t,0) 



(G6) 



p(t, x) 



OO 

E 



e rj p y ' 



(t,x) 



(Gl) 



where tl c = (M/R 3 ) 1 / 2 . 

At order 0(e B r] 1 ) the fluid equations are the same as 
Eqs. (|G2f) with (0,2) — > (1,1) and the driving accelera- 
tion replaced by 



a) = —\v 



and similar equations for P, $, and v. These expansions 
are plugged into the fluid equations Q3. l|l and solved at 
each order in e and 77. We set e = eg and ignore electric- 
type tidal interactions. 15 

At order 0{e 1 B rp) for n < 2, the results are the same as 
in Sec. IIIIBl with the relabelling (n) — > (n, 0). At order 
0( £ e 7 ? 1 ) the velocity perturbation is unconstrained and 
chosen to be uniform rotation, in ' 1 ) = Si x x. The 
density and other fluid variables are unchanged at this 
order: = P^ 1 ) = S^ 1 ) = 0. 

At order 0(e%rj 2 ), we find the usual decrease in central 
density due to rotation. This can be calculated explic- 
itly using the methods of Sec. IIII Bl The perturbation 
equations at this order become 

dt 

and 



^ x B]i 



,(1,0) 



V]v. 



(0,1) 



lij k 3$ % 1 



(G7) 



where 



^ijk — ^ y^^jk^i $ijQk ^aj^a^ik ^abj^-ick^ac^b) ■ 

(G8) 

Since the angle average of n ■ a^ 1,1 ) vanishes, there is no 
change in central density at this order. [If we consider 
electric-type tidal interactions, there should also be no 
change in central density at order 0(e 1 £ rf L ). This follows 
from the fact that it is impossible to construct a scalar 
that is linear in both and Qj.] 

At order 0(e B rj 2 ), the driving acceleration becomes 



+ V^°'%f 2) ] = , (G2a) M = _ [t>(0>1) . ^ a,!, _ [v(M) . v]v 



(0,1) 



(G9) 



,(1,0) 



(0,2) 



,(0,2) 



V\V. 



(1,0) 



b (0 < 2) x B]i 



Jo, 2 



dt p(°.°) 



1 vy (^(0,2) 1 P ' vy <f,(o,o) _ „(o,2) At each order, the Lagrangian perturbation of the 



,(0,0) 



star is approximately given by Q n,m ^ + w? 
fn,m)^ go ^ e velocity perturbations scale like t^™'" 1 - 1 cx 



*(n,m) 



(G2b) 



where 



a (0,2) = _ [t>(0l l) , VK (0,1) = _ {n _ a , )Q . _ (G3) 

The change in central density is then determined by con- 
verting to an equation for the Lagrangian displacement 
of the fundamental radial mode fSec. lIII C|l . The angular 
integral of n - a^ 0,2 ) is (8/3)7rrf2 2 , and the resulting inner 
product with £0 is 



(£o,a 



(0,2) 



0.4163Mi? 2 f7 2 , 



(G4) 



where we have computed the following radial integral for 
a r = 2 polytrope using the techniques in appendix iFl 



u 2 £(u) sinudu w 14.43 



(G5) 



15 Combined electric- type/magnetic- type tidal couplings cannot 
change the central density up to order 0(a 7 Q 2 ). Any changes to 
the central density up to this order must have the form 8p c /p c ~ 
EijBij h £ abc £ ad B bd n c h £ ab B bc Q a n.c, where "&" means "plus 
terms of the form". Parity considerations force all these terms 
to vanish: Under a parity transformation (xj — » — ay), Sp c /p c — > 



Spc/po, while £i- 



£ij , Bi 



-Bij, £l{ > rii, and €ijk 



—f-ijk- All three of the above terms are parity odd while the 
change in central density is parity even. 



{n,m) 1^2^ density perturbations scale like 

p(n,m) ^ _p(o,o)yi i f(™, m )^ rp determine if the change in 
central density vanishes at a given order, one can com- 
pute the angle average of n ■ a^ n,rn \ Since each of the 
terms in n ■ a/ 1 ' 2 ) is proportional to an odd power of rij, 
their angular integrals vanish and one finds that there 
is no change in central density at order 0(e B r/ 2 ). One 
can also argue that at order 0(e B r] 2 ) the change in cen- 
tral density must be proportional to BijSliQj and must 
vanish because it is parity odd. (If electric-type tidal 
interactions were included, a change in central density 
proportional to Sijfliflj would be allowed.) 

Following the same procedure at order 0(e B rj 1 ), one 
finds the total acceleration 

af < 1} = -[vW . v] v f'°) - [„d.°) • Vh (M) (G10) 
-[v^ ■ Vh (0 ' 1} - [v^ ■ Vh (2 < 0) + [vW x B]t . 

In this case, the change in central density must be propor- 
tional to E a h c B a dBbdQc, which does not obviously vanish 
from parity arguments. A change in central density at 
this order is therefore possible. If electric-type tidal inter- 
actions are considered, a central density change at order 
0{e 2 : , q l ) proportional to e a b c £ad£bd^c is also allowed. 
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At order 0(e|?7 2 ) the central density must be propor- 
tional to n 2 BijBij or BijBik^lj^lk- Similarly, for electric- 
type tidal fields the change in central density at order 
0(e|?7 2 ) must be proportional to Vl 2 £ij£ij or c', ,'',/. <! ,<!/. . 
These terms have even parity and need not vanish. 

To summarize, we have shown by simple arguments 



that in rotating stars the change in central density van- 
ishes at orders O^ry 1 ), 0{e%rr x ), 0(e^°), O^rf), 
0{e 1 £ rr x ), 0{e 1 B rr x ), andO(e^?7 2 ). Determining the change 
in central density at other orders would require more ex- 
plicit calculations. 
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